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Parameter  Set  Estimation  of  Time  Varying 
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By 
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Stephen  Yurkovich,  Adviser 


Parameter  set  estimation  (PSE),  a  class  of  system  identification  schemes  which 
aim  at  characterizing  the  uncertainty  in  the  identification  experiment,  will  play  a  vital 
role  in  robust  identification  for  control.  An  important  step  in  current  research  along 
these  lines  is  development  of  PSE  algorithms  for  systems  which  are  time  varying  in 
nature;  this  is  particularly  true  if  the  identified  model  set  is  to  be  used  in  an  adaptive 
setting. 

In  this  dissertation,  the  Optimum  Volume  Ellipsoid  (OVE)  algorithm  for  pa¬ 
rameter  set  estimation  of  time-invariant  systems  is  extended  to  time-varying  sys¬ 
tems.  Building  on  this  development  of  the  OVE  algorithm  for  Time- Varying  systems 
(OVETV),  two  algorithms  are  also  presented  for  reducing  the  computational  com¬ 
plexity  of  the  optimal  time  update  equations.  These  algorithms,  scalar  addition  and 
scalar  multiplication,  reduce  the  computational  complexity  of  the  time  update  equa- 
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tions  by  constraining  the  new  ellipsoid  to  be  parameterized  by  the  previous  ellipsoid 
and  a  single  new  parameter.  Several  examples  are  presented  in  detail. 

Following  this,  extensions  to  the  OVE  and  OVETV  algorithms  are  presented.  The 
algorithms  are  extended  to  handle  multiple-input,  single-output  (MISO)  systems.  It 
is  shown  how  knowledge  of  dependencies  in  the  parameter  variations  can  be  exploited 
to  reduce  the  number  of  computations  in  the  resulting  algorithm.  A  “square-root” 
implementation  of  the  OVE  algorithm  is  developed  which  has  improved  numerical  sta¬ 
bility  properties.  We  show  how  the  scope  of  OVE-ISP,  an  OVE-based  input  synthesis 
procedure,  can  be  extended  to  handle  systems  with  known  transportation  lag. 

Lastly,  it  is  shown  how  the  OVE  and  OVETV  algorithms  can  be  utilized  for  fault 
detection  and  isolation  (FBI).  Two  methods  are  suggested  for  detecting  faults  in 
dynamical  systems.  The  first  method  relies  on  a  consistency  check  which  is  integral 
to  the  OVE  and  OVETV  algorithms,  while  the  second  method  utilizes  an  ellipsoid 
intersection  test  to  detect  a  fault.  Two  “recovery”  strategies  are  also  presented  which 
allow  the  OVE  and  OVETV  algorithms  to  continue  to  track  the  parameters  after  a 
fault  is  detected. 

During  the  course  of  this  research,  we  derived  some  extremely  useful  routines  using 
symbolic  computations  for  linear  time-varying  systems;  these  appear  in  an  appendix 
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CHAPTER  I 


Introduction 


1.1  Motivation 

Throughout  the  centuries,  humans  have  desired  to  control  their  environment.  Within 
the  last  150  years,  a  large  body  of  theory  has  been  developed  for  designing  algorithms 
which  are  capable  of  controlling  elements  of  this  environment  autonomously.  Most 
of  these  algorithms  depend  either  explicitly  of  implicitly  on  a  mathematical  model 
of  the  system  to  be  controlled.  There  are  two  fundamental  methods  for  developing 
models  which  describe  the  dynamic  behavior  of  a  system. 

The  first  method  is  an  analytical  approach  whereby  these  mathematical  equations 
are  developed  from  the  basic  laws  of  physics,  chemistry,  biology,  etc.  One  drawback 
with  this  approach  is  the  fact  the  system  may  be  too  complex  for  the  derivation  of  a 
model  to  be  practical  or  even  possible  using  this  approach. 

In  contrast,  the  second  method  for  deriving  a  model,  system  identification,  is  “rel¬ 
atively’  easy  to  use.  System  identification  is  an  empirical  approach  whereby  experi¬ 
ments  are  performed  on  the  system,  and  parameter  values  of  a  model  are  estimated 
based  on  the  observed  data.  This  method  is  not  without  its  drawbacks.  First  of  all,  it 
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Parameter  (Set) 
Estimate 


Figure  1:  System  identification  scheme 

may  not  be  valid  for  all  inputs  and  operating  points.  Secondly,  in  general,  it  does  not 
provide  as  much  physical  insight  into  the  problem  as  the  analytical  approach.  Most 
complex  systems  usually  require  some  combination  of  these  approaches. 

In  spite  of  these  drawbacks,  system  identification  is  a  very  powerful  tool  and  the 
subject  of  this  dissertation.  A  typical  block  diagram  of  a  system  identification  scheme 
is  shown  in  Figure  1.1,  where  the  plant  is  the  system  to  be  modeled,  u  is  the  system 
input,  and  y  is  the  system  output.  The  estimator  attempts  to  find  parameters  (or 
a  set  of  parameters)  which  describe  the  behavior  of  the  system.  Unfortunately,  this 
process  is  complicated  by  the  fact  that  this  must  be  done  in  the  presence  of  noise  and 
uncertainty,  v. 
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Not  surprisingly,  the  characterization  of  v  dictates  the  nature  of  the  algorithm 
used  to  estimate  the  parameters.  Most  traditional  parameter  estimation  algorithms 
assume  that  v  is  a  random  variable  which  satisfies  certain  statistical  properties.  Given 
these  assumptions,  the  algorithms  seek  to  find  a  “best”  estimate  of  the  parameters 
that  describe  the  system.  A  problem  with  this  approach  is  the  fact  that  this  “best” 
estimate  will  probably  not  match  the  “true”  parameters  exactly.  Consequently,  there 
is  no  guarantee  that  a  control  algorithm  designed  for  this  “best”  estimate  will  perform 
satisfactorily  on  the  actual  system. 

Alternatively,  it  can  be  assumed  that  v  is  unknown-but-bounded.  This  results  in  a 
parameter  set  which  is  theoretically  guaranteed  to  contain  the  true  parameter.  If  the 
control  algorithm  is  designed  for  the  entire  parameter  set,  some  sense  of  robustness 
has  been  gained. 

Conceptually,  parameter  set  estimation  is  a  simple  idea.  Consider  the  system 

yk  =  buk  +  Vk,  (1-1) 

where 

-l<Vk<l,  (1.2) 

and  where  b  is  the  parameter  we  seek  to  estimate.  Solving  (1.1)  for  Vk  and  substituting 
into  (1.2)  results  in  —1  <  yk  —  buk  <  1  or 

yk-l  <buk<yk  +  l-  (1.3) 

An  experiment  was  performed  on  the  system  from  which  the  data  in  Table  1  was 


measured. 
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Table  1:  Measured  data 


k 

Uk 

Vk 

1 

0.5 

1.5 

2 

O 

r— 1 

1.5 

Substituting  the  measured  values  at  A;  =  1  into  (1.3),  we  find  that  the  parameter  b 
must  lie  within  a  certain  feasible  range,  or  set,  of  parameters;  that  is,  b  must  satisfy  b  € 

[1.5] .  Likewise,  at  A:  =  2,  6  must  also  satisfy  b  6  [0.5, 2.5].  The  parameter  b  is  constant 
(time-invariant),  and  therefore,  must  satisfy  both  these  conditions.  Therefore,  the 
actual  value  of  b  must  be  contained  in  the  “feasible  parameter  set”  [1,2.5]  where 

[1.2.5]  =  [1,5]  n  [0.5, 2.5]. 

Unfortunately,  as  this  concept  is  extended  to  higher  dimensions,  the  feasible  pa¬ 
rameter  set  is  quite  complicated  and  difficult  to  track  computationally.  Because  of 
their  computational  efficiency  and  ease  in  mathematical  expression,  ellipsoids  are 
often  used  to  overbound  the  feasible  parameter  set.  However,  if  the  estimated  pa¬ 
rameter  set  is  to  be  used  for  robust  control,  there  is  a  tradeoff  between  the  volume 
of  the  ellipsoid  and  the  achievable  performance  level.  Consequently,  we  would  like  to 
find  the  smallest  volume  ellipsoid  which  is  guaranteed  to  contain  the  true  parameter 
values. 

This  was  the  motivation  for  the  Optimal  Volume  Ellipsoid  (OVE)  algorithm  which 
was  developed  in  [1]  for  single-input,  single-output  linear  time-invariant  systems.  This 
algorithm  is  computationally  efficient  and  is  guaranteed  to  contain  the  true  plant. 
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While  the  OVE  algorithm  was  developed  for  linear  time-invariant  systems,  there 
is  a  large  number  of  systems  where  the  underlying  model  to  be  identified  is  linear 
time- varying.  The  identification  of  time- varying  systems  is  an  integral  part  of  many 
adaptive  control  and  adaptive  signal  processing  applications.  Consequently,  the  pri¬ 
mary  goal  of  this  dissertation  will  be  to  extend  the  OVE  parameter  set  estimation 
algorithm  for  identification  of  time-varying  systems. 


1.2  Review  of  the  Literature 

While  the  work  in  this  dissertation  was  motivated  in  particular  by  [1],  there  has  been 
a  large  amount  of  research  in  this  area.  Most  parameter  set  estimation  algorithms 
use  either  ellipsoids  or  boxes  to  bound  the  feasible  parameter  set;  one  exception  is 
found  in  [2]  where  the  exact  feasible  parameter  set  is  calculated. 

The  most  well  known  methods  for  bounding  the  feasible  parameter  set  are  the 
Optimal  Bounding  Ellipsoid  (OBE)  algorithms.  Fogel  and  Huang  [3]  were  the  first 
to  apply  overbounding  ellipsoids  to  the  parameter  set  estimation  (PSE)  problem.  In 
their  algorithm,  the  center  of  the  bounding  ellipsoid  can  be  viewed  as  a  particular 
weighted  least  squares  estimate.  They  considered  two  optimization  strategies.  One 
strategy  sought  to  minimize  the  volume  of  the  ellipsoid  while  the  second  strategy 
sought  to  minimize  the  lengths  of  the  ellipsoid  semi-axes. 

A  third  optimization  strategy  with  no  interpretable  geometrical  significance  was 
used  by  Dasgupta  and  Huang  in  [4].  In  [5],  it  was  shown  how  these  optimization 
strategies  can  be  combined  with  different  weighting  strategies  to  include  a  larger 


6 


CHAPTER  I.  Introduction 


family  of  algorithms  that  they  call  unified  OBE  (UOBE).  Included  in  this  family  is 
the  set-membership  weighted  recursive  least  squares  (SM-WRLS)  algorithm. 

These  algorithms  are  recursive  in  nature.  They  also  have  the  capability  of  ignor¬ 
ing  redundant  data.  However,  even  when  they  use  the  minimum  volume  optimization 
strategy,  they  do  not  produce  the  ellipsoid  with  the  minimum  volume.  This  conser¬ 
vatism  is  undesirable  for  robust  control  design.  Although  the  volume  convergence  of 
the  Fogel  and  Huang  algorithm  was  improved  by  Belforte  and  Bona  [6],  it  still  did 
not  produce  an  ellipsoid  with  the  minimum  volume  because  the  center  estimate  was 
constrained  to  be  a  weighted  least  square  estimate. 

In  [1],  Cheung,  Yurkovich,  and  Passino  eliminated  this  constraint  which  led  to  the 
Optimal  Volume  Ellipsoid  (OVE)  algorithm,  shown  to  have  desirable  convergence 
properties  (proof)  and  geometrical  interpretation.  Two  independently  developed  al¬ 
gorithms,  which  have  been  shown  to  be  mathematically  equivalent  to  the  OVE  algo¬ 
rithm,  are  the  the  modified  OBE  (MOBE)  algorithm  of  [7]  and  the  EPC  algorithm 
of  [8].  The  difference  between  the  algorithms  lies  in  convergence  proofs,  geometrical 
interpretation,  and  implementation. 

Several  extensions  to  the  OVE  algorithm  have  already  been  considered.  In  [9], 
a  modified-OVE  (MOVE)  algorithm  was  presented  which  allows  for  time-varying 
bounds  on  v  and  simplifies  the  parameter  update  calculations.  The  OVE  algorithm 
was  extended  for  interconnected  systems  in  [10].  In  [11],  an  OVE-based  input  syn¬ 
thesis  procedure  was  presented  which  seeks  to  choose  the  system  input  such  that  a 
rapid  decrease  in  ellipsoid  volume  occurs. 
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We  are  particularly  interested  in  how  PSE  can  be  extended  to  time-varying  sys¬ 
tems.  As  in  the  time-invariant  case,  most  adaptive  PSE  algorithms  have  used  over¬ 
bounding  ellipsoids.  However,  two  algorithms  using  polytopes  can  be  found  in  [12].  In 
[13],  three  ellipsoid-type  algorithms  for  tracking  time  varying  systems  are  discussed: 
scalar  bound  inflation,  fixed-memory  bounding,  and  bound  incrementing.  Most  adap¬ 
tive  PSE  algorithms  fit  into  one  of  these  three  categories. 

The  simplest  of  these  is  scalar  bound  inflation.  In  this  algorithm,  the  ellipsoid  is 
expanded  uniformly  in  all  direction  before  processing  the  new  measurement.  While 
computationally  simple,  one  difficulty  with  this  approach  has  been  choosing  the  ex¬ 
pansion  factor  such  that  the  set  estimate  is  guaranteed  to  contain  the  true  plant. 
Similar  to  the  forgetting  factor  in  the  standard  weighted  recursive  least  squares  al¬ 
gorithm,  the  expansion  factor  is  often  chosen  heuristically,  usually  to  be  a  constant 
[5]. 

In  [14],  Rao  and  Huang  show  how  the  Dasgupta- Huang  OBE  algorithm  can  be 
modified  by  artificially  increasing  the  noise  bound  to  guarantee  consistency  if  an 
update  takes  place.  However  if  no  update  takes  place,  this  consistency  can  still  be 
lost.  They  also  have  developed  a  procedure  to  “recover”  when  tracking  is  lost,  for 
example  due  to  an  unexpected  large  parameter  jump. 

With  fixed  memory  bounding  the  data  is  windowed,  that  is,  only  the  last  L  obser¬ 
vations  are  used  to  find  the  plant  estimate.  Norton  and  Mo  [13]  found  this  method 
awkward  to  implement  with  no  advantage  over  scalar  bound  inflation  or  bound  incre¬ 
menting.  However,  in  [5]  “back- rotation”,  an  efficient  method  for  implementing  the 
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window,  is  applied  to  SM-WRLS.  With  this  strategy  two  additional  windowing  tech¬ 
niques  were  also  considered.  In  the  first,  instead  of  a  rectangular  window,  tapered 
windowing  is  used  where  the  window  tapers  slowly  to  zero.  In  the  second,  selec¬ 
tive  forgetting  is  used  where  only  previously  heavily  weighted  data  sets  are  removed. 
These  algorithms  offer  no  guarantee  of  consistency. 

The  last  strategy  that  [13]  considered  was  the  bound  incrementing  method.  This 
method  is  based  on  an  explicit  parameter  variation  model  in  which  the  possible  pa¬ 
rameter  deviation  during  each  time  step  is  assumed  to  be  ellipsoidally  bounded.  They 
used  a  result  given  in  [15]  which  finds  the  minimum  volume  ellipsoid  that  bounds  the 
sum  of  two  vectors  which  are  also  ellipsoidally  bounded. 

A  significant  motivation  for  extending  PSE  to  time-varying  systems  is  adaptive 
robust  control.  As  defined  in  [16],  the  robust  control  problem  is  the  problem  of  ana¬ 
lyzing  and  designing  accurate  control  systems  given  plants  which  contain  significant 
uncertainty.  While  this  is  still  an  open  research  area,  there  has  been  some  work  in 
this  direction. 

In  [17],  Kosut  examined  some  of  the  issues  involved  in  using  PSE  for  adaptive 
control.  A  finite  horizon  controller  was  designed  for  a  finite  impulse  response  (FIR) 
model  whose  parameters  were  ellipsoidally  bounded  in  [18].  In  [19],  a  controller  was 
designed  for  systems  whose  output  matrix  was  ellipsoidally  bounded.  Bound-based 
worst-cased  self-tuning  controllers  were  designed  in  [20]. 

Some  very  promising  results  for  robust  control  can  be  found  in  [21].  In  this  work, 
techniques  are  developed  for  extracting  modeling  uncertainty  information,  vital  for 
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robust  control  design,  from  the  results  of  PSE  algorithms.  A  worst-case  analysis 
was  used  to  address  the  characterization  of  parametric  uncertainty  as  additive  and 
coprime  factor  perturbations  to  a  nominal  plant  model.  Techniques  for  identification 
of  reduced-order  perturbation  weightings  were  also  given. 


1.3  Dissertation  Organization 

This  dissertation  is  organized  into  four  main  chapters.  In  Chapter  II,  the  OVE 
algorithm  is  extended  for  identification  of  time-varying  systems.  After  illustrating 
how  the  equations  arise  for  an  optimal  time  update  when  there  are  two  parameters 
to  estimate,  the  results  of  [15]  are  combined  with  the  optimum  volume  measurement 
update  equations  of  [1]  to  develop  an  Optimal  Volume  Ellipsoid  algorithm  for  Time- 
Varying  systems  (OVETV). 

Building  on  the  development  of  the  OVETV,  two  algorithms  are  presented  for 
reducing  the  computational  complexity  of  the  optimal  time  update  equations.  These 
algorithms  reduce  the  computational  complexity  of  the  time  update  equations  by 
constraining  the  new  ellipsoid  to  be  parameterized  by  the  previous  ellipsoid  and  a 
single  new  parameter.  One  of  these  algorithms  uses  the  same  parameterization  as 
the  scalar  bound  inflation  method  of  [13].  However,  the  algorithms  developed  in  this 
chapter  are  optimized  for  volume  and  combined  with  the  optimal  volume  measurement 
update  equations  of  [1]. 

Three  examples,  which  illustrate  the  main  feature  of  the  OVETV  algorithm,  are 
considered  in  Chapter  III.  The  first  is  a  simple  first-order  example  which  compares 
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the  OVETV  algorithm  with  the  scalar  addition  and  scalar  multiplication  approaches. 
In  the  second  example,  identification  of  a  linear  time- varying  circuit  is  used  to  demon¬ 
strate  how  this  approach  can  be  applied  to  sampled-data  systems  with  quantization 
noise.  In  the  final  example,  the  OVETV  algorithm  is  applied  to  actual  data  obtained 
from  a  crude  oil  distillation  column. 

In  Chapter  IV,  four  extensions  to  the  OVETV  algorithm  are  presented.  In  many 
real  systems,  the  output  of  the  system  is  effected  by  more  than  one  input.  We  discuss 
how  the  OVE  and  OVETV  algorithms  can  be  used  to  identify  multiple- input,  single¬ 
output  (MISO)  systems.  Results  are  demonstrated  on  the  distillation  column  example 
of  Chapter  III. 

Often,  the  time- varying  parameters  that  we  are  trying  to  identify  are  dependent  on 
a  smaller  set  of  physical  time- varying  parameters.  If  we  can  exploit  this  dependence, 
we  may  be  able  to  simplify  the  problem  and  reduce  the  number  of  computations 
required  to  implement  the  algorithm.  Adaptations  of  the  OVETV  algorithm  are 
developed  to  utilize  knowledge  of  these  dependencies. 

When  applying  the  OVE  and  OVETV  algorithms  in  practice,  numerical  round-off 
errors  can  cause  the  algorithms  to  become  numerically  intractable.  To  solve  this  prob¬ 
lem,  we  apply  an  approach  that  has  been  successfully  used  in  stochastic  estimation 
schemes.  A  “square-root”  algorithm  is  developed  for  implementing  OVE.  Further¬ 
more,  it  is  shown  how  the  square-root  approach  can  be  extended  to  time-varying 
systems. 

The  type  of  system  input  for  identification  experiments  has  a  very  significant 
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effect  on  the  “quality”  of  estimates  produced  from  an  identification  scheme.  The 
final  extension  to  the  OVETV  algorithm  which  we  will  consider  in  Chapter  IV  is 
based  on  the  OVE-ISP  procedure  of  [11].  We  show  how  the  scope  of  OVE-ISP  can 
be  extended  to  handle  systems  with  known  transportation  lag  such  as  the  distillation 
column  in  Chapter  III.  We  also  demonstrate  its  application  to  time-varying  systems 
and  an  important  “stabilizing”  property  that  other  synthesis  strategies  do  not  appear 
to  possess. 

Fault  detection  and  isolation  (EDI)  is  concerned  with  the  detection  and  identifi¬ 
cation  of  failures  in  complex  dynamical  systems.  In  Chapter  V,  we  discuss  how  the 
OVE  and  OVETV  algorithms  can  be  utilized  for  FDI.  Two  methods  are  suggested 
for  detecting  faults  in  dynamical  systems.  The  first  method  relies  on  a  consistency 
check  which  is  integral  to  the  OVE  and  OVETV  algorithms,  while  the  second  method 
utilizes  an  ellipsoid  intersection  test  to  detect  a  fault. 

Two  “recovery”  strategies  are  presented  which  allow  the  OVE  and  OVETV  algo¬ 
rithms  to  continue  to  track  the  parameters  after  a  fault  is  detected.  The  first  strategy 
simply  resets  the  current  ellipsoid  to  an  ellipsoid  which  is  “large  enough”  to  guar¬ 
antee  that  the  “true”  parameter  is  contained  in  the  parameter  set  immediately  after 
a  fault  is  detected.  An  alternative  strategy  based  on  projections  is  also  introduced. 
While  this  algorithm  is  not  guaranteed  to  recapture  the  ellipsoid,  it  usually  produces 
ellipsoids  with  smaller  volumes. 

An  integrated  approach  is  also  suggested  for  combining  these  recovery  strategies. 
Finally,  the  algorithms  of  Chapter  V  are  illustrated  using  the  linear  time-vaxying 
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circuit  example  of  Chapter  III. 

Finally,  during  the  course  of  this  research,  investigation  into  time- varying  linear 
systems  was  imperative.  In  so  doing,  we  derived  some  extremely  useful  routines  using 
symbolic  computations  for  linear  time-varying  (LTV)  systems.  To  our  knowledge, 
such  routines  have  not  appeared  elsewhere  in  the  open  literature;  we  have  gathered 
this  material  in  a  stand-alone  unit,  given  in  Appendix  B. 


CHAPTER  II 


OVE  for  Time- Varying  Systems 

2.1  Overview 

Identification  of  time- varying  systems  is  essential  for  many  adaptive  control  and  adap¬ 
tive  signal  processing  applications.  In  particular,  adaptive  tracking  algorithms  are 
needed  for  (i)  complex  systems  which  admit  a  linear,  time- varying  parameterization, 
(ii)  failure  detection  and  identification  (FDI)  systems,  and  (iii)  gain  scheduling  con¬ 
trol  systems  where  supervisory  controllers  must  track  changing  system  characteristics. 
While  most  adaptive  tracking  algorithms  assume  that  the  perturbations  to  the  system 
satisfy  certain  statistical  properties,  in  many  real  world  applications  it  is  often  more 
natural  to  assume  that  these  perturbations  are  unknown-but-bounded. 

Traditional  system  identification  techniques  abound  for  such  applications,  where 
a  point  in  the  parameter  set  is  identified  to  model  the  system.  A  different  philosophy, 
discussed  in  [22],  [17],  and  elsewhere,  is  to  seek  to  identify  a  set  of  models  which 
are  consistent  with  the  incoming  data  (referred  as  set  membership  identification  or 
parameter  set  estimation,  or  PSE).  A  major  motivation  for  adopting  this  philosophy 
for  the  identification  portion  of  the  overall  control  design  process  is  for  use  in  robust 
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control  design.  While  many  adaptive  aad  auto-tuning  techniques  have  been  developed 
based  on  point  estimation  schemes  over  the  last  20  years,  it  has  only  been  recently 
that  connections  have  been  made  to  the  robust  control  design  problem  for  ellipsoidally 
based  PSE  techniques;  see,  for  example,  [21]  and  [19].  It  is  for  this  reason,  within  this 
promising  line  of  research  for  robust  identification  and  control,  that  we  focus  here 
on  PSE  algorithms,  and  limit  our  discussions  and  comparisons  with  point  estimation 
schemes.  Indeed,  there  exists  a  large  body  of  literature,  which  we  cannot  address 
here,  on  point  estimation  schemes,  many  of  which  are  useful  in  practice  (see,  e.g., 
[23]). 

Generally  speaking,  PSE  algorithms  are  classified  as  robust  identification  tech¬ 
niques,  because  they  attempt  to  characterize  uncertainty  by  identifying  the  set  of 
feasible  parameters,  given  the  data  from  the  identification  experiment.  That  is,  PSE 
seeks  to  identify  a  set  of  parameters  which  are  feasible  with  the  measured  data  and 
the  bounds  on  the  perturbations.  In  general,  this  feasible  set  is  an  irregular  convex 
set.  Because  of  their  computational  efficiency  and  ease  in  mathematical  expression, 
ellipsoids  are  used  to  overbound  the  feasible  set.  When  PSE  is  used  in  robust  control 
or  robust  adaptive  control,  there  is  a  tradeoff  between  set  size  and  system  perfor¬ 
mance.  Consequently,  one  seeks  the  smallest  volume  ellipsoid  which  is  guaranteed  to 
contain  the  “true”  plant  parameters. 

In  [13],  the  authors  considered  the  bound  incrementing  and  scalar  bound  infla¬ 
tion  methods  for  PSE  of  time- varying  systems.  In  [1],  the  Optimal  Volume  Ellip¬ 
soid  (OVE)  algorithm  was  developed  for  a  time-invariant  single-input,  single-output 
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(SISO)  auto-regressive  with  exogenous  input  (ARX)  model,  and  extended  in  [10]  for 
interconnected  systems.  The  bounded  incrementing  method  of  [13]  uses  an  optimum 
volume  algorithm  found  in  [15]  for  updating  the  time  equations. 

In  this  chapter,  after  illustrating  how  the  equations  arise  for  the  optimal  time 
update  equations  when  there  are  two  parameters  to  estimate,  we  will  combine  the 
results  of  [15]  with  the  optimum  volume  measurement  update  equations  of  [1]  to 
develop  an  Optimal  Volume  Ellipsoid  algorithm  for  Time- Varying  systems  (OVETV). 
The  OVETV  enjoys  many  of  the  favorable  characteristics  of  the  OVE,  which  in  fact 
are  shared  with  other  mathematically  equivalent  techniques  such  as  extensions  of  the 
modified  OBE  algorithm  of  [7],  and  the  ellipsoid  with  parallel  cuts  (EPC)  algorithm 
of  [8]. 

Building  on  the  development  of  the  OVETV,  we  present  two  algorithms  for  re¬ 
ducing  the  computational  complexity  of  the  optimal  time  update  equations,  thereby 
making  real-time  implementation  feasible  for  many  applications.  These  algorithms 
reduce  the  computational  complexity  of  the  time  update  equations  by  constraining 
the  new  ellipsoid  to  be  parameterized  by  the  previous  ellipsoid  and  a  single  new 
parameter.  One  of  these  algorithms  uses  the  same  parameterization  as  the  scalar 
bound  inflation  method  of  [13].  However,  the  algorithms  developed  in  this  chapter 
are  optimized  for  volume  and  combined  with  the  optimal  volume  measurement  update 
equations  of  [1]. 
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2.2  Problem  Statement 

Consider  the  time- varying  SISO  ARX  model 


Vk 

=  Ok<l>k  +  Vk, 

(2.1) 

6kJr\ 

=  6k  +  Wk. 

(2.2) 

where  6k  =  an  • 

'  '  ^nk  ^Ik 

iT 

bmk  is  the  parameter  vector,  yk  is 

the  system 

output,  Uk  is  the  system  input,  and  <i>k  =  [  — j/fc-i  •  •  •  —yk-n  Uk-i  •  •  Uk-m 
is  the  regression  vector.  Here  I  is  the  number  of  time  steps  delay  in  the  system  where 
0  <  /  <  m.  Let  r  =  n  +  m  —  l-Il.  The  system  disturbance,  Ufc,  and  the  parameter 
disturbance  vector,  Wk  E  are  assumed  to  satisfy  the  following  bounds 


2k<vk<  ')k, 

(2.3) 

wlRl^Wk  <  1, 

(2.4) 

where  Xk  <  ^k  Xki  ^k  G  are  known  at  each  time  k.  The  vector, 

Wk,  represents  the  change  in  the  parameter  vector  at  each  time  k.  The  matrix,  Rk, 
is  symmetric,  positive  definite  and  represents  an  ellipsoidal  bound  on  the  possible 
parameter  vector  deviation  during  each  time  step. 

Equations  (2.1)  and  (2.2)  and  the  bounds  in  (2.3)  and  (2.4)  are  used  to  define 
the  adaptive  parameter  set  estimation  problem.  Let  Hk-i  C  be  the  set  such  that 
all  Ok-i  €  Hk-i  are  feasible  parameters  of  the  plant  which  are  consistent  with  past 
Vk,  <t>ki  If.,  lki  ^k-  ^  Let  Fk  C  3^’’  be  the  set  such  that  all  6k  E  Fk  are  feasible 

^In  this  work,  the  notation  A  C  B  means  that  A  is  a  subset  of  B,  and  therefore,  A  may  be  equal 
to  B. 
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parameter  estimates  of  the  plant  which  are  consistent  with  the  measurements  at  time 
k.  That  is,  from  (2.1)  and  (2.3), 

Fk  =  {Ok  €  31?’'  :  7fc  <  2/fc  -  Ol<i>k  <  7jt},  (2.5) 


where  Fk  is  the  region  between  two  parallel  hyperplanes. 

Let  Gk  C  3?’’  be  the  set  such  that  all  Ok  €  Gk  are  feasible  parameter  estimates  of 
the  plant  which  are  consistent  with  Ok-i  G  iffc-i-  That  is,  from  (2.2)  and  (2.4), 

Gk  =  [Ok  e  r  :  {Ok  -  Ok-ifRkUOk  -  Ok-i)  <  l;0fc-i  €  (2.6) 

Before  incorporating  a  new  measurement,  the  set,  Gk,  tells  us  how  much  the  parame¬ 
ter  set  estimate  at  time  {k  —  1),  Hk-i,  must  expand  to  guaxantee  consistency  at  time 
k.  For  Hk  to  be  consistent,  it  must  satisfy  (2.5)  and  (2.6);  therefore,  Hk  is  given  by: 

Hk  =  GknFk.  (2.7) 

The  problem  of  parameter  set  estimation  for  time-varying  systems  is  to  find  Hk  ex¬ 
plicitly  in  the  parameter  space  where  Hk  is  defined  recursively  by  (2.5),  (2.6),  and 
(2.7).  In  general,  however,  computing  (2.7)  is  an  extremely  diflScult  problem  because 
the  feasible  parameter  space,  Hk,  is  an  irregular  convex  set.  Consequently,  ellipsoids 
will  be  used  to  overbound  Hk  because  of  their  ease  in  mathematical  expression  and 
computational  efficiency. 

Define  the  ellipsoid  Ek-i  as 


Ek-i  =  [Ok-i  :  (0;b-i  -  h-ifPFM&k-i  -  h-i)  <  l-,0k-i  e  r}  (2.8) 
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where  9k-i  G  3^’’  is  the  center  of  the  ellipsoid,  and  Pk-i  G  is  a  symmetric, 

A 

positive  definite  matrix.  Similarly,  define  Ek  as 

Ek  =  {Ok  :  (Ok  -  hfpp\0k  -  h)  <  1;  G  r }  (2.9) 

A  ^ 

where  9k  G  3?’'  is  the  center  of  the  ellipsoid  and  Pk  G  is  a  symmetric,  positive 
definite  matrix. 

Choose  Ek-i  such  that  Ek-\  D  Hk-i-  By  replacing  Hk-i  with  Ek-i  in  (2.6),  we 

A 

obtain  a  new  set,  Gk-,  given  by 

Gk  =  {9k  G  r  :  {9k  -  9k-ifRklA^k  -  9k-i)  <  1;  9k-i  G  Ek-i}.  (2.10) 

A  AAA 

In  general,  Gk  is  not  an  ellipsoid.  Consequently,  we  choose  Ek  such  that  Ek  D  Gk, 
i.e.,  we  wish  to  bound  the  set,  Gk  with  an  overbounding  ellipsoid.  In  particular,  we 
would  like  to  choose  the  Ek  with  minimum  volume,  that  is 

Ek  =  arg£;  min{vol(£)  :  E  D  Gk}-  (2-11) 

In  addition  to  Gk,  at  time  k  the  ellipsoid,  Ek,  is  also  constrained  by  Fk  as  defined  in 
(2.5).  Clearly,  we  wish  to  choose  Ek  such  that  Ek  D  Ek  Ci  Fk,  i.e.,  we  wish  to  bound 

A 

all  consistent  parameter  estimates  with  an  overbounding  ellipsoid.  As  with  Ek,  we 
seek  the  Ek  with  minimum  volume,  that  is 

Ek  =  axgE  niin{vol(£l)  :  E  D  EkC  Fk}-  (2.12) 

Equations  (2.5),  (2.10),  (2.11),  and  (2.12)  provide  a  recursive  algorithm  for  finding 
an  ellipsoid  Ek  such  that  Ek  D  Hk.  However,  this  algorithm  assumes  that  the  solutions 
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to  (2.11)  and  (2.12)  can  be  calculated.  The  solution  to  (2.12),  the  measurement 
update  equations,  will  be  solved  using  a  modified  version  of  the  OVE  algorithm 
found  in  [9].  The  solution  to  (2.11),  the  time  update  equations,  will  be  developed  for 
r  =  2,  i.e.,  when  there  are  two  parameters  to  estimate.  A  result  by  Chernous’ko  will 
be  used  to  generalize  this  to  the  full  r-dimensional  case  [15]. 


2.3  Measurement  Update  Equations 

In  [1],  Cheung,  Yurkovich,  and  Passino  developed  the  OVE  algorithm,  with  conver¬ 
gence  proofs,  for  linear  time  invariant  systems,  that  is,  when  Wk  =  0.  This  result  was 
extended  for  interconnected  systems  in  [10].  The  OVE  algorithm  solves  the  following 
optimization  problem: 

Ek  =  arg£;  min{vol(E)  :  E  3  Ek-i  0  Fk}.  (2.13) 

The  most  significant  difference  between  the  OVE  algorithm  and  the  seminal  work  of 
[3]  is  that  unlike  the  OVE  algorithm,  the  OBE  algorithm  constrains  the  center  of  the 
new  ellipsoid,  to  satisfy  a  “modified”  recursive  least  squares  estimate.  It  is,  in 
fact,  true  that  the  OVE,  the  modified  OBE  (MOBE)  algorithm  of  [7],  and  the  EPC 
algorithm  of  [8]  are  mathematically  equivalent.  The  difference  between  the  algorithms 
lies  in  convergence  proofs,  geometrical  interpretation,  and  implementation. 

The  OVE  algorithm  in  [1]  assumes  that  7^  =  —7^  =  7.  In  [9],  Gassman  and 
Yurkovich  discuss  the  modified  OVE  (MOVE)  algorithm  which  uses  the  more  general 
bound  in  (2.3).  The  optimization  in  (2.13)  is  identical  to  (2.12)  with  Ek  replaced  by 
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Ek-i .  Consequently,  we  will  be  able  to  use  the  MOVE  algorithm  by  simply  replacing 

^  A  A 

dk-i  and  Pk-i  with  9k  and  Pk,  respectively. 

With  these  substitutions,  the  algorithm  is  given  as  follows: 

1.  Set 


If  Qjj.  >  1  or  Ofe  < 
algorithm  stops. 

2.  Set  tk  =  OfcOfe.  If  Cfc  <  —  then  no  measurement  update  is  necessary,  that  is 

~dk  =  0k  A  =  Pk.  (2.16) 


(Vk-Glh-ik 
.  ,yk-H<f>k-ik 

au  =  mini - — 1 ) 


(2.14) 


(2.15) 


-1,  then  the  observed  data  is  inconsistent  with  Ek  and  the 


3.  Set  fik  =  ^2'~- 

4.  If  l/ifcl  >  p,  then 


bk 

o  1  +  Cfc 

=  2rpk  + 

Pk 

(2.17) 

Tk 

bk  -  s[gn{pk)\/bk  -  4(r  +  1)(1  +  rek) 

2(r  +  1) 

(2.18) 

(Tk 

=  Tk{Tk  )  +  1 

Pk 

(2.19) 

bk 

(Tk 

l_Ik 

(2.20) 
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5.  If  \pk\  <  P,  then  =  0  and 


a  =  max(|a^.|,  \ak\) 

CTk  =  ro? 

.  r(l  —  a^) 


6.  Update  6k  and  Pk 


6k 

Pk 


6  k  + 


TkPk<f>k 

i<j>lPk<Pky/^ 


^kPk  +  {<^k—  6k) 


Pk<t>k<j>kPk 

<PkPk<f>k 


(2.21) 

(2.22) 

(2.23) 


(2.24) 

(2.25) 


To  arrive  at  the  equations,  an  affine  transform  can  be  used  to  transform  the 
ellipsoid  Ek  to  a  unit  ball  centered  at  the  origin.  In  the  transformed  coordinate 
system,  a*,  and  are  the  coordinates  along  <f)k  of  the  two  hyperplanes  (defined  by 
Fk)  which  are  orthogonal  to  (j>k.  The  consistency  check  in  step  1  is  used  to  verify  that 
the  intersection  of  Ek  and  Fk  is  not  empty.  If  this  intersection  is  empty,  then  either 
Eo  did  not  contain  the  “true”  parameter  or  the  assumptions  made  in  (2.1)  through 
(2.4)  were  invalid.  Recovery  strategies  have  been  developed  to  handle  the  situation 
where  an  inconsistency  does  occur. 

The  check  in  step  2  is  used  to  determine  whether  the  current  data  record  contains 
enough  information  to  reduce  the  ellipsoid  volume.  In  steps  4  and  5,  p  is  chosen  to 
be  a  very  small  value  and  is  used  to  determine  when  w  —  a^.  In  the  transformed 
coordinate  system,  t^,  <f)k,  and  8k  are  used  to  characterize  the  optimum  volume  ellip¬ 
soid  Ek-  The  parameter  Tk  is  the  center  coordinate  of  Ek  along  (j>k,  (Jk  is  the  squared 
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length  of  the  semi-axis  of  Ek  along  (f)k,  and  6k  is  the  squared  length  of  the  semi-axes 
of  Ek  orthogonal  to  (f)k.  In  step  6,  the  center  Ok  and  orientation  matrix  Pk  of  the 
optimum  volume  ellipsoid  Ek  are  given  in  the  original  coordinate  system.  Finally,  it 
should  be  noted  that  a  similar  algorithm  for  set  membership  state  estimation  can  be 
found  in  [24]. 

2.4  Time  Update  Equations 

Initially,  we  will  constrain  Rk  to  be  where  I  is  the  identity  matrix,  and  Cfc  bounds 
||wjk||2.  Doing  so,  equation  (2.10)  can  be  rewritten  as 

Gk  =  {Ok  e  r  :  {Ok  -  0k-if{0k  -  Ok-x)  <  Cl-i;  Ok-i  e  Ek-i).  (2.26) 

Later,  this  constraint  will  be  relaxed. 

Clearly  for  Cfc-i  =  0  where  k  =  I ...  N ,  Ek  =  Gk  =  Ek-i,  that  is,  this  algorithm 
reduces  to  the  MOVE  algorithm.  Note  that,  while  Eo  cannot  be  chosen  to  be  3J’',  it 
can  be  chosen  arbitrarily  large  so  that  the  “true”  parameter  vector,  Oq^  is  guaranteed 
to  be  in  Eq. 

A 

To  solve  (2.11),  it  is  important  that  we  understand  the  nature  of  Gk-  From 
(2.26),  we  see  that  Gk  is  a  union  of  hyperspheres  with  radius  (k  whose  centers  are 
contained  in  the  ellipsoid,  Ek-i-  The  hyperspheres  of  greatest  interest  are  those  on 

A 

the  hypersurface  of  the  ellipsoid  because  they  define  the  hypersurface  of  Gk-  See 
Figure  2  for  an  example  of  the  2-D  case. 

Let  z7  be  the  unit  vector  normal  to  the  hypersurface  of  Ek  that  points  away  from 
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a 


Figure  2:  Union  of  Ek-i  with  spheres  on  the  surface  of  Ek-i 

the  interior  of  Ek-  Clearly,  the  distance  along  iy  between  the  hypersurface  of  Ek  and 
the  hypersurface  of  Gk  is  always  Cfc-i-  In  fact,  this  defines  the  hypersurface  of  Gk- 
The  surfaces  of  Ek-\  and  Ek  are  defined  by 

=  {Si-i  :  -  h-ifPtlAh-i  -  h-i)  =  1;  «t-i  €  »’),  (2.27) 

Si  =  :  (h  -  hfPi'iek  -  h)  =  1;  9i  €  S'}.  (2.28) 

Let  the  distance  between  two  parameter  vectors,  Oi  and  Oj,  be  given  by 

d{9i,ej)=\\ei-9j\U.  (2.29) 


Next  we  need  the  following  lemma. 
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A 

Lemma  1  The  minimum  distance  between  the  hypersurface  of  Ek  and  the  hypersur¬ 
face  of  Ek-x  must  he  Cit-i;  that  is, 

Cfe_i  =  min{c?(4_i,^fc) :  Ok-x  €  Sk-x,Ok  €  Sk}  (2.30) 


A  A 

Proof;  If  the  minimum  distance  is  less  than  Cfc-i,  then  Ek  f)  Gk-  This  is  a  contradic¬ 
tion  because  Ek  must  bound  Gk-  If  the  minimum  distance  is  greater  than  (k-x,  then 
there  must  exist  some  ellipsoid  E  D  Gk  where  vol(£?)  <  vol{Ek)-  This  is  a  contradic- 
tion  because  Ek  must  be  the  minimum  volume  ellipse  that  contains  Gk-  Therefore, 
(2.30)  must  hold.  □ 

Using  Lemma  1,  equation  (2.11)  can  be  rewritten  as 

Ek  =  arg£;min{vol(£')  :  min(d(0fc_i, 0fc))  =  (k-x',0k-i  G  Sk-x,0k  G  (2.31) 

That  is,  we  want  to  find  the  minimum  volume  ellipsoid,  Ek,  such  that  the  minimum 

A 

distance  between  the  surface  of  Ek-x  ^-nd  Ek  is  equal  to  Cfc-i-  This  suggests  a  two  step 
optimization  problem.  First,  solve  (2.30)  analytically  for  minimum  values  Ok-x  and 
Ok-  That  is,  find  Ok-y  and  6^  as  a  function  of  Ek-x  and  Ek  such  that  d{9k_i,0k)  =  Cfc-i- 
Secondly,  substitute  9l_x  and  91  into  (2.31)  and  solve  for  the  minimum  volume  Ek  as 
a  function  of  Ek-x- 

Using  the  singular  value  decomposition  (SVD),  the  positive  definite  matrix,  Pk-x, 
used  to  define  the  ellipsoid  in  (2.8)  can  be  factored  as 

Pk-x  =  Vk-xAk-xVlx 


(2.32) 
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where  Vk-\  €  is  orthogonal  and  Afc_i  G  is  diagonal.  The  axes  of  the  ellipsoid 
are  aligned  with  the  columns  of  Vk-\ .  The  semi-axes  lengths  are  given  by  the  diagonal 
elements  of  [25].  Note  also  that  vol(A?fe_i)  oc  \Jdei[Pk-\)  —  y^det(Afc_i).  Due  to 

A  A 

the  symmetry  of  Ek-i  and  Gk,  the  minimum  volume  Ek  must  have  the  same  center 
and  axes  as  Ek,  that  is, 


k  =  h-i,  (2.33) 

Pk  =  Vk-ikkVl^.  (2.34) 

Consequently,  the  optimization  problem  becomes  one  of  finding  A*  as  a  function  of 
Afc_i  such  that  (2.30)  is  satisfied  and  ^det(A*;)  is  minimized. 

With  this  said,  the  general  r  dimensional  problem  will  be  discussed  later.  The 
two  dimensional  problem,  r  =  2,  is  solved  next. 


2.4.1  Two-Dimensional  Case 

Although  most  systems  that  we  would  like  to  identify  have  more  than  two  parameters 
to  estimate,  the  two-dimensional  problem  is  important  because  it  is  easy  to  visualize 
and  it  provides  insight  into  the  r-dimensional  problem.  In  the  two-dimensional  case, 
“ellipsoids”  become  “ellipses”  and  “volume”  is  actually  “area”. 

As  discussed  in  the  previous  section,  we  will  begin  by  solving  (2.30)  for  and 
0\.  By  (2.33)  and  (2.34),  it  is  clear  that  we  can  consider  ellipses  centered  at  the  origin 
with  axes  aligned  along  the  coordinate  axes  without  any  loss  of  generality.  Define  the 
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surface  of  ellipses,  Ei  and  E2,  as 

2  2 

-^1  =  {a:i,2/i  :  ^  +  ^  =  e  3ft},  (2.35) 

®i 

2  2 

<5'2  =  {a;25  2/2 :  ^  ^  =  l;a:2,j/2  €  3ft},  (2.36) 

«2 

where  ai,  fei,  02,  and  62  are  positive  constants.  The  distance  between  Si  and  S21 
(2.29),  can  be  rewritten  as 

d{xi,x2,yuy2)  =  \/{xi  -  X2y  +  (yi  -i/2)^-  (2.37) 

Now  we  can  prove  the  following  lemma.  Note,  that  because  of  symmetry,  we  can  find 
the  solution  in  the  upper  right  quadrant  without  any  loss  of  generality. 

Lemma  2  The  parameters  which  minimize  (2.37),  the  distance  between  Si  and  S2, 


must  satisfy  one  of  four  solutions: 


•  Minimum  distance  along  y-axis: 

x\  =  0  yl  =  hi,  (2.38) 

x*2  =  0  yl  =  62,  (2-39) 

with 

d{xl,xl,yl,y2)  =  \b2-bi\,  (2.4O) 

•  Minimum  distance  along  x-axis: 

x*i  =  ai  yl  =  0,  (2.4I) 

X*2  =  02  y*i=  0,  (3.42) 
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with 

d{xl,xl,yi, 

•  Ellipses  Intersect: 

x*i  =Pa 
X*2  =P4 

with 

d(x’,X2. 


;)  =  |02  -  Oil,  (&PJ 

</,*  =  Ps,  (i-U) 

vl  =  Pi.  (US) 

i:,P2')  =  0,  (t46) 


•  Minimum  distance  between  axes: 


.  «1\/^ 

KVp: 

(U7) 

"  (a?  -  bl)^ 

(<■!  -  4?)VS’ 

*  «2V^ 

(•ivK 

(US) 

-  (ai  -  62)^ 

(4  -  6i)v^’ 

with 


where 


d{x\,xl,yl,yl)  = 

{b\al  -  a\bl){b\  af} 

(8.49) 

(62  —  02)(6i  —  a{) 

P\  =  a\hl-h\al, 

(2.50) 

P2  =  —  2620201  +  O2O1  —  61‘*a2  +  2h\a\a\  —  aV^a\, 

(2.51) 

P3  =  h\h\  —  6j62  —  2h\h\a\  +  6201  +  2h\h\a\  —  a^b^, 

(2.52) 

aia2Jbl  -  bl 

P4  =  ", 

^/Pi 

(2.53) 

6162^0?  -  al 

Ps  =  ^  • 

(2.54) 
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Proof:  Because  d{xi,X2,yi,y2)  is  nonnegative,  minimizing  cP  is  equivalent  to  mini¬ 
mizing  d.  Using  the  Lagrange  multiplier  method,  the  constraints  posed  by  the  ellipses, 
(2.35)  and  (2.36),  can  be  augmented  to  d?  to  form 

f(xi,x2,yi,y2,^i,^2)  =  {xi  -  X2y  +  {yi  —  y2Y  +  (2.55) 

+ylal  -  albl)  -f  X2{xlbl  +  ylal  ~ 

For  f{xl,X2,yl,y2,Xiy\l)  to  be  a  minimum,  the  differential,  df{xl,xl,yi^y2,  XI,  X^) 
must  be  0.  Taking  the  differential  of  (2.55),  we  get  the  following  set  of  equations  that 
must  be  solved 


2{xl  -  xl) Xl{2xlbl)  =  0, 

(2.56) 

Wi  -  yl)  +  =  0. 

(2.57) 

2{xl-x;)-i-X*{2x*2bl)  =  0, 

(2.58) 

m  -  y2)  +  Ki‘^yW2)  =  0, 

(2.59) 

+  =  0) 

(2.60) 

xfbl  +  yfal  -  albl  =  0. 

(2.61) 

The  solution  of  (2.56)  through  (2.61)  yields  Lemma  2.  □ 

The  second  step  of  the  optimization  process  requires  substitution  of  xj,  y^,  x^, 
and  2/2  into  (2.31).  Doing  so  and  solving  for  the  minimum  volume  ellipse,  E2,  results 
in  the  following  lemma. 


Lemma  3  The  minimum  volume  ellipse,  E2,  whose  surface,  S2,  is  at  least  C  >  0 
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from  the  surface  Si  of  ellipse  Ei  is  parameterized  by 

a  =  I  ^  ai  =  hi 
\  V&  otherwise 

[  ^\/P8  +  (a?  -  bl)^  ai  >  bi 

^2  =  <  6i  +  C  ai  =  bi 

[  7^\/Ps  +  (^1  -  ai)\/^  «1  < 

where  pi  through  ps  are  defined  in  Appendix  A. 


(2.62) 

(2.63) 


Proof:  As  mentioned  above,  we  begin  by  substituting  xj,  yl,  x%,  and  y^  into  (2.31). 
However,  Lemma  2  gives  four  possible  solutions  which  depend  on  the  values  of  oi,  6i, 
02,  and  62-  Clearly,  we  are  not  interested  in  the  3rd  case,  that  is,  when  the  ellipses 
intersect.  The  case  that  we  are  most  interested  in  is  case  4,  that  is,  when  the  minimum 
distance  is  between  the  axes.  However,  case  4  is  not  a  solution  when  oi  =  61. 

We  begin  by  considering  the  case  oi  =  and  assume  that  the  minimum  distance 
is  along  the  j/-axis,  i.e.,  case  1.  We  wish  to  minimize  the  volume  of  E2,  which  is 
equivalent  to  minimizing 

y  =  0262.  (2.64) 

If  the  minimum  distance  is  along  the  y-axis  then  62  —  61  =  C  +  C-  Also 

implied  by  (2.40)  is  that  02  —  01  >  C-  Consequently,  02  =  oi  +  C  if  (2.64)  is  to  be 
minimized.  Clearly,  if  we  assume  that  the  minimum  distance  was  along  the  x-axis, 
we  would  get  the  same  result. 

Now  consider  the  case  when  Oi  7^  bi  and  assume  that  the  minimum  distance  is 
between  the  axes.  The  constraint  posed  by  (2.49)  can  be  appended  to  (2.64)  giving 


K(a2.  i2,  A)  =  oA  +  A((6l  -  0=  -  4?  +  -  44)  -  4(4  -  4)i4  -  a?))  (2.65) 
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to  be  minimized.  For  Va{al,bl,  \*)  to  be  a  minimum,  the  differential,  dVa{al^b2,X*), 
must  be  0.  Taking  the  differential  of  (2.65),  we  get  the  following  set  of  equations  that 
must  be  solved 

b*2  +  2X*{alalbf  -  26^af  +  bja^bf  +  bjalal  -  b^a^ 

-C^a^al  +  C^albl)  =  0  (2.66) 

a*2  +  2A*(6;6^af  -  2aj6f  +  alb^al^  +  alb^bl  -  a^b^ 

-CVi  +  C%al)  =  0  (2.67) 

(6f  -  -bl  +  al){bla*2^  -  albf)  -  (^bf  -  a*'^){bl  -  af)  =  0  (2.68) 

The  solution  of  (2.66),  (2.67),  and  (2.68)  yields  Lemma  3.  □ 

We  can  now  prove  the  following  theorem. 

Theorem  1  For  r  =  2,  let  Pk-i  he  factored  as  Pk-\  =  where  Vk-\  G 

is  orthogonal  and 

A).-.=  j,  .  (s  e9) 

Then,  the  minimum  volume  ellipse,  Ek,  which  contains  the  set,  Gk,  where  Ek  is 
defined  by  (2.9)  and  Gk  is  defined  by  (2.26),  is  given  by 

Pk  =  Vk-ikkVl,  (2.70) 

where 

A.=  [f  “].  (H-71) 

and  02  and  h^  are  given  by  Lemma  3  with  (  replaced  by  (k-i  • 
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Proof:  The  proof  follows  directly  from  equations  (2.33)  and  (2.34)  and  Lemmajs  1-3. 

□ 


Example  of  Minimum  Volume  Ellipse 


In  this  section  we  will  apply  Theorem  1  to  a  particular  Ek-i  and  (k-i-  Let  (^k-i  =  0.5 
and  Ek-i  be  given  by  (2.8)  with 


Pk-i  = 


25  -5 
-5  5  J  ’ 

6  ■ 

7  ■ 


(2.72) 

(2.73) 


The  SVD  of  Pk-i  is  given  by 


-0.9732  -0.2298 
0.2298  -0.9732  ’ 


(2.74) 


Afc-1 


26.1803  0 

0  3.8197 


(2.75) 


From  (2.69),  we  get  ci  and  bi  which  are  used  by  Lemma  3  to  solve  for  02  and  62- 
Substituting  02  and  62  into  (2.71)  results  in 


L  = 


32.6125  0 

0  6.1300 


Using  (2.33)  and  (2.34),  Ek  is  given  by  (2.9)  with 


(2.76) 


Pk 


31.2145  -5.9217 
-5.9217  7.5279  ’ 


6 

7 


(2.77) 


(2.78) 
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a 


Figure  3:  Minimum  volume  ellipse,  Ek,  that  bounds  Gk 

In  Figure  3  we  see  that  Ek  does  bound  Gk-  Note  that  although  Ek  has  the  smallest 
volume  of  any  ellipse  that  contains  Gjt,  it  is  still  larger  than  Gk  except  for  the  case 
oi  =  bi.  However,  except  for  the  cases  where  oi  '>  bi  or  ai  <C  bi  this  difference  in 
volume  is  small. 

2.4.2  r-Dimensional  Case 

To  generalize  this  to  r-dimensions,  we  will  use  a  result  by  Chernous’ko.  In  [15], 
Chernous’ko  finds  the  minimum  volume  ellipsoid  which  bounds  the  sum  of  two  vectors 
which  are  also  ellipsoidally  bounded.  This  is  exactly  the  problem  posed  by  (2.11). 
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The  solution  is  given  as  follows: 


1.  Solve  the  generalized  eigenvalue  problem 

for  each  Aj,  j  €  [l,r]. 

2.  Solve 

’'1  r 

Y' _ f _ = _ 1 _ 

+P  PiP+ 1) 

for  the  unique  p  >  0. 

3.  Then,  the  solution  to  (2.11)  is  given  by 

^k  = 

Pk  =  {p  +  i)Pk-l  +  {p  ^  +  1)-Ra:-1 


(2.79) 


(2.80) 


(2.81) 

(2.82) 


To  arrive  at  these  equations,  a  nonsingular  matrix  can  be  found  which  transforms 
the  ellipsoid  Ek-i  to  a  unit  ball  and  the  ellipsoid  Gk  to  an  ellipsoid  whose  semi¬ 
axes  are  aligned  with  the  coordinate  axes.  After  the  transformation,  the  semi-axes 
length  of  the  ellipsoid  Gk  are  given  by  the  square  root  of  the  generalized  eigenvalues 
of  (2.79).  The  optimization  in  (2.11)  is  then  solved  in  the  transformed  coordinate 
system.  This  optimization  requires  the  solution  of  (2.80).  A  proof  that  (2.80)  has 
only  one  positive  root  can  be  found  in  [15].  After  transforming  back  to  the  original 
coordinate  system,  the  solution  to  (2.11)  is  given  by  (2.81)  and  (2.82).  Note  that  Rk 
is  no  longer  constrained  to  be  C|/. 
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2.5  OVETV  Algorithm  Summary 

The  steps  of  the  Optimal  Volume  Ellipsoid  for  Time- Varying  systems  (OVETV)  al¬ 
gorithm  are  summarized  below. 

1.  Choose  the  initial  ellipsoid,  Eq,  “large  enough”  such  that  the  initial  “true” 
parameter  vector,  is  in  Eq. 

2.  Find  Ek  such  that 

jEfc  =  arg£;min{vol(jE)  :  £  D  Gfc}.  (2.83) 

3.  Find  Ek  such  that 

Ek  =  arg£:  min{vol(£')  :  E  D  EkC  Fk}.  (2.84) 

4.  Repeat  steps  two  and  three  for  each  new  measurement. 

2.6  Alternative  Strategies 

The  optimal  time  update  equations  require  the  solution  of  an  rth  order  generalized 
eigenvalue  problem  (2.79)  and  an  (r  -1- 1)  order  polynomial  (2.80).  Computationally, 
this  may  not  be  feasible  for  certain  real-time  applications,  especially  when  r  is  large. 
Consequently,  we  will  consider  two  alternative  algorithms  which  represent  different 
degrees  of  tradeoff  between  computational  complexity  and  resulting  volume.  These 
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algorithms  also  use  the  MOVE  algorithm  to  solve  (2.12);  thus,  we  will  only  discuss 
the  approximating  solutions  to  (2.11).  These  algorithms  reduce  the  computational 
complexity  by  parameterizing  the  new  ellipsoid  with  one  parameter  instead  of  r  pa¬ 
rameters.  Again,  we  will  begin  with  the  two-dimensional  case. 


2.6.1  Two-Dimensional  Case 

If  we  consider  the  2-dimensional  system  described  in  Section  2.4.1,  one  possible  pa¬ 
rameterization  is  given  by 


7]ai 

(2.85) 

rjbi. 

(2.86) 

The  parameterization  in  (2.85)  and  (2.86),  scalar  multiplication,  is  appealing  because 
it  does  not  require  finding  all  the  singular  values  of  Pk-i.  Note  that  this  parameteri¬ 
zation  was  refered  to  as  scalar  bound  inflation  in  [13].  Unfortunately,  it  may  yield  a 
very  conservative  bounding  ellipsoid.  For  Ek  to  contain  Gk,  the  following  conditions 
on  the  semi-axes  must  hold: 


>  «i  +  C 

(2.87) 

>  &1  +  C. 

(2.88) 

Substituting  (2.86)  into  (2.88)  and  solving  for  rj,  we  get 


(2.89) 
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Substituting  (2.89)  into  (2.85)  yields 

«2  ^  +  Ci~’  (2.90) 

h 

Clearly,  for  oi  '>  bi,  02  can  become  very  large. 

For  an  example,  see  Figure  4.  Here,  we  set  Oi  =  9,  61  =  1,  and  vary  (  from  0  to  1.9. 
This  plot  contains  02  and  62  which  are  the  optimal  values  calculated  using  Theorem  1. 
Also  plotted  are  the  lower  bounds,  (2.88)  and  (2.90),  for  the  parameterization  in  (2.85) 
and  (2.86).  Clearly,  for  C  >  0  this  parameterization  will  produce  an  ellipse  which  is 
much  larger  than  optimum.  Despite  this  conservativeness,  we  will  state  the  following 
lemma  without  proof.  It’s  implications  will  be  discussed  in  the  next  section. 

Lemma  4  Let  a2  =  ait],  62  =  bit]  where  rj  >  0.  Let  the  surfaces  of  ellipses,  Ei  and 
E2,  be  defined  by  (2.35)  and  (2.36),  then  the  value  oft]  which  minimizes  the  volume 
of  E2  such  that  the  surface,  S2,  is  at  least  C  >  0  from  Si  is 

V  =  l  +  -,  (^M) 

X 

where  x  =  min(ai,  61). 

Examining  Figure  4,  it  appears  that  the  change  in  02  versus  (  and  the  change  in  62 
versus  C  are  similar.  This  suggests  parameterizing  the  new  ellipse  by  scalar  addition, 
that  is 

«2  =  ai-\-T}  (2.92) 

&2  =  bi-\-T],  (2.93) 

which  leads  to  the  following  lemmas. 
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Figure  4:  Comparison  of  ellipse  parameters  between  optimum  and  scalar  multiplica¬ 
tion 
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Lemma  5  Let  a2  =  ai-\-T],  62  =  + »?  where  rj  >0.  Let  the  surfaces  of  ellipses,  E\ 

and  E2,  be  defined  by  (2.35)  and  (2.36).  Then  the  minimum  distance  between  Si  and 


S2  is  given  by 


d=J2y^- 


2biai  +  bip  +  aip 


(61  +  ai){bi  +2t]  +  ai) 

Proof:  This  lemma  follows  from  Lemma  2.  Substituting  (2.92)  and  (2.93)  into  (2.49) 
yields  (2.69).  A  more  direct  proof  will  be  used  to  extend  this  method  to  r  dimensions. 


□ 


Lemma  6  Let  02  =  «i  +  »?;  ^2  =  where  y  >  0.  Let  the  surfaces  of  ellipses, 

El  and  E2,  be  defined  by  (2.35)  and  (2.36).  Then  the  value  ofy  which  minimizes  the 
volume  of  E2  such  that  the  surface,  S2,  is  at  least  C  >  0  from  Si  is 

»  =  PT  +  P6  -  1^,  (S.95) 

where 


Po 

=  cii  bi 

(2.96) 

Pi 

=  276t  +  366?Gi  +  2aj6^  +  366ia^  +  2a^ 

(2.97) 

P2 

=  i^3(-16pgC"+PiC2 
Po  ^ 

-  646?af ) 

(2.98) 

P3 

-albl 

Spl 

(2.99) 

P4 

= 

(2.100) 

P5 

2po 

(2.101) 

Pe 

,/  2  1  64 

=  Ws*’"  27'^ 

(2.102) 
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Proof:  By  Lemma  1,  for  E2  to  have  minimum  volume,  the  minimum  distance,  d,  in 
(2.94)  must  be  equal  to  (.  Solving  (2.94)  for  t}  yields  (2.95).  □ 

We  can  now  prove  the  following  theorem. 


Theorem  2  For  r  =  2,  let  Pk-i  be  factored  as  Pk-i  =  where  Vk-i  € 

is  orthogonal  and 

A/.-.  =  jj  .  (S.IO4) 

Let  the  semi-axes  of  Ek  be  parameterized  as  02  =  ai  +  rj  and  62  =  61  +  »/.  Then,  the 
minimum  volume  ellipse,  Ek)  which  contains  the  set  Gk  where  Ek  is  defined  by  (2.9) 
and  Gk  is  defined  by  (2.26),  is  given  by 


Pk  =  Vk-iAkVl, 


(2.105) 


where 

h  =  I  m 

and  7}  is  given  by  Lemma  6  with  (  replaced  by  (k-i  • 

Proof:  The  proof  follows  directly  from  equations  (2.33)  and  (2.34)  and  Lemmas  5 
and  6.  □ 

In  order  to  compare  the  values  of  Theorem  2  with  the  optimum  values  of  Theo¬ 
rem  1,  again  we  set  oi  =  9  and  61  =  1.  This  time  (  is  varied  from  0  to  8.9.  Shown  in 
Figure  5  are  the  plots  of  02  and  62  from  Theorem  1  and  Theorem  2.  Clearly,  this  is 
a  much  closer  approximation  than  the  one  given  by  (2.85)  and  (2.86).  The  volumes 
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Figure  5:  Comparison  of  ellipse  parameters  between  optimum  and  scalar  addition 


were  calculated  for  the  parameters  shown  in  Figure  5.  A  ratio  of  the  volume  using 
scalar  addition  versus  the  optimum  volume  is  shown  in  Figure  6.  For  this  exam¬ 
ple,  the  volume  using  scalar  addition  is  at  most  approximately  7%  higher  than  the 
optimum  volume. 


Theorem  2  presents  a  reasonable  alternative  to  Theorem  1.  In  the  next  section, 
we  will  extend  it  as  well  as  scalar  multiplication  to  r  dimensions. 
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Figure  6:  Ratio  of  scalar  addition  volume  versus  optimum  volume 


2.6.2  r-Dimensional  Case 

Scalar  Addition 

In  the  extension  to  r-dimensions,  we  see  that  Ek  is  constrained  to  be  an  ellipsoid  with 
the  same  center  and  orientation  as  Ek-\  and  with  semi-axis  lengths  which  are  found 
by  adding  a  scalar,  rj,  to  the  semi-axis  lengths  of  Ek-\.  Rk  is  again  constrained  to  be 
Cfc/;  therefore,  Gk  is  given  by  (2.26). 

As  in  the  2- dimensional  case,  we  need  to  parameterize  the  minimum  distance 
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between  the  following  ellipsoids.  Define  the  surface  of  ellipsoids,  Ei  and  E2,  as 

=  (2-107) 


.=1 


^2  =  {5Gr:X] 


r 

X; 


+  vY 

where  each  ai,  i  €  [l,r]  and  t]  are  positive  scalars. 


=  1} 


(2.108) 


Lemma  7  Let  the  surfaces  of  ellipsoids,  E\  and  E2,  be  defined  by  (2.107)  and  (2.108). 
Then,  the  minimum  distance  between  Si  and  S2  is  given  by 


2aa  +  ay  +  ay 
(a  +  a)(a  +  2?;  +  a) 


(2.109) 


where 


a  =  max  a,-  a  =  mm  a,-. 

j6[l,r]  ^  -  j6[l,r] 


(2.110) 


Proof:  Because  d{x,x)  is  nonnegative,  minimizing  cP  is  equivalent  to  minimizing  d. 
Using  the  Lagrange  multiplier  method,  the  constraints  posed  by  the  ellipsoids,  (2.107) 
and  (2.108),  can  be  augmented  to  (P  to  form 


T  ^2 

X; 


fix,  X,  h,  A,)  =  J:(x,  -  +  A,(S  -  1)  +  HY,  7-^  -  1).  (2.111) 

i=l  i=l  “i  t=l  "r  'Ij 

For  /(a:*,x*,  Aj,  Aj)  to  be  a  minimum,  the  differential,  df(x*,  £*,  Aj,  A2)  must  be  0. 
Taking  the  differential  of  (2.111),  we  get  the  following  set  of  equations  that  must  be 
solved 


X 


2{x]-x])  +  2\l-^  =  0  i€[l,r] 


-  2(x*j  -  X*)  +  2A* 


X: 


(«i  +  vY 


=  0  j  e  [l,r 


(2.112) 

(2.113) 
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=  0 


i=i 


=  0 


(2.114) 

(2.115) 


^  («i  +  v)^ 

The  solution  of  (2.112)  through  (2.115)  yields  Lemma  7.  □ 

This  lemma  holds  because  if  the  semi-axes  of  ellipsoid,  Bi,  are  incremented  by  a 
scalar  t],  then  the  minimum  distance  between  the  surface  of  Ei  and  the  surface  of 
the  resulting  ellipsoid,  E2,  must  lie  along  the  two-dimensional  cross-section  defined 
by  the  maximum  and  minimum  length  semi-axes. 

We  now  state  the  following  theorem: 

Theorem  3  Let  Ek-i  and  Ek  be  defined  as  in  (2.8)  and  (2.9).  Let  Gk  be  given  by 
(2.26).  If  the  SVD  of  Pk-\  is  given  by  Pk-\  =  T4_iAfc_i\4T.i  where  Vk-i  is  orthogonal 
and  Ak-i  =  diag{\i . . .  A^),  then  the  value  of  y  which  optimizes 


rj  =  arg,  min{t;o/(E)  :EDGk;P  =  Vk-^A]!!^  -H 


is  given  by 


where 


2  aa 


a  =  max  \  \i  a  =  min  \  \i, 

je[i,ri  V  “  icn.riV 


(2.116) 


(2.117) 


(2.118) 


and 


Po  =  a  + a 

Pi  =  27a^  H-  36a^a  -|-  2a^a^  -t-  36a^  +  2a* 


(2.119) 

(2.120) 
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P2  = 

^±±^Z{-liaCU  +th<U  -64a3#) 

Po  ^ 

(2.121) 

P3  = 

—orgr 

spi 

(2.122) 

P4  = 

-\Ck-iPo 

(2.123) 

P5  = 

mCk-i 

2po 

(2.124) 

P6  = 

,/  2  1  64 

y-?'  -  5P4 + 27P3  - 

(2.125) 

pr  = 

./  2  1  _^64 

y-jPs  -  jw  +  + 

(2.126) 

Proof:  The  proof  follows  directly  from  equations  (2.33)  and  (2.34)  and  Lemmas  6 
and  7.  □ 

Scalar  Multiplication 

In  the  examples  we  have  considered,  Theorem  3  provides  an  acceptable  alternative  to 
the  optimal  algorithm.  Notice  that  the  (r  +  1)  order  polynomial  in  (2.80)  has  been 
reduced  to  the  same  computational  complexity  as  the  2nd  order  problem.  However, 
the  algorithm  still  requires  the  solution  of  an  rth  order  SVD. 

With  this  in  mind,  we  reconsider  the  scalar  multiplication  strategy  discussed  in 
section  2.6.1.  In  this  algorithm,  Ek  is  constrained  to  be  an  ellipsoid  with  the  same 
center  and  orientation  as  Ek-i;  the  semi-axis  lengths  are  found  by  multiplying  the 
semi-axes’s  of  Ek-i  with  a  scalar,  rj. 

Theorem  4  Let  Ek-i  and  Ek  be  defined  as  in  (2.8)  and  (2.9).  Let  Gk  be  given  by 
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(2.10).  Then  the  value  of  rj  which  optimizes 

r]  =  arg^  mm{vol{Ek)  :  Ek  D  Gk, Ok  =  4-i,  A  =  V^h-i}  (2.127) 

is  given  by 

i]  =  l  +  \fi  (2.128) 

where  A  is  the  maximum  generalized  eigenvalue,  Xj,  which  satisfies 

Rk-iXj  =  XjPk-iXj.  (2.129) 

Proof:  From  [26],  we  know  that  E  D  Gk  il  and  only  if 


SeSO)  >  (2-130) 

where  S{6)  defines  the  support  function  of  the  respective  sets.  We  also  know  that 

=  +  (2.131) 

with  the  support  functions  given  by 

SE,iO)  =  O^h  +  y/o^,  (2.132) 

Se,.^{0)  =  e'^ek-i  +  \/oTpi^,  (2.133) 

Su„_M  =  y/0TRk-iO.  (2.134) 

Substituting  (2.131)-(2.134)  into  (2.130)  results  in 

e'^Ok  +  \/9^PkO  >  O^h-I  +  ^O^Pk-iO  +  y/o'^Rk-iO.  (2.135) 


By  substituting  6k  =  A-i  and  A  =  fpPk-\  into  (2.135),  it  is  easy  to  show  that 


(7/  -  1)"^ A-1^  >  O'^Rk-iO. 


(2.136) 
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From  [27],  we  know  that  there  exists  a  matrix  transformation,  0  =  T6,  such  that 


{ti  -  ifFe  > 


(2.137) 


where  A  is  a  diagonal  matrix  containing  the  generalized  eigenvalues  of  (2.129).  For 
(2.137)  to  hold,  it  is  necessary  that  (r/  -  1)^  >  A  where  A  is  the  maximum  generalized 
eigenvalue  of  (2.129);  thus  rj  >  1  +  \/A.  The  volume  of  Ek  is  proportional  to  ^et(A:) 
which  is  minimized  if  the  equality  holds;  therefore,  7/  =  1  +  V\.  □ 

Note  that  this  algorithm  does  not  require  the  solution  of  an  (r  +  1)  order  polyno¬ 
mial.  Furthermore,  the  rth  order  generalized  eigenvalue  problem  of  (2.79)  is  replaced 
with  finding  the  largest  eigenvalue.  There  exist  iterative  techniques  which  can  be 
used  to  find  this  value  quickly.  We  could  also  use  the  relationship,  ||Pfc"\i2fc-i||  >  A 
[28],  to  replace  77  in  (2.128)  with  77  =  1  -b  where  ||  •  ||  is  any  matrix  norm; 

note  that  ^  >  77.  Also,  note  that  Pk-i  is  given  recursively  by  [29] 


4  "  Vfc  6k’(j>lPk<f>k 


(2.138) 


where  P^  ^  =  ^Pk\- 

Theorem  4  always  gives  more  conservative  results  than  the  OVETV .  For  the  case 
where  Rk  =  CkP  Theorem  4  often  gives  more  conservative  results  than  Theorem  3. 
However,  the  reduced  number  of  computations  may  make  it  the  most  acceptable  for 
real  time  applications. 
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2.7  Summary 

In  this  chapter,  we  developed  the  optimal  volume  ellipsoid  algorithm  for  parameter 
set  estimation  of  time-varying  systems.  We  also  developed  two  alternative  strategies 
that  offer  various  degrees  of  trade  off  between  computation  complexity  and  ellipsoid 
volume.  Given  our  assumptions  in  (2.1)-(2.4),  all  three  algorithms  are  guaranteed 
to  contain  the  “true”  parameter  in  the  parameter  set.  While  OVETV  gives  the 
optimum  volume  ellipsoid,  the  scalar  addition  and  scalar  multiplication  algorithms 
give  optimum  volume  ellipsoids  assuming  additional  constraints  on  the  time  update 
ellipsoid. 


CHAPTER  III 


Examples 


3.1  Overview 

In  his  chapter,  we  consider  three  examples.  The  first  is  a  simple  first-order  exam¬ 
ple  which  serves  to  illustrate  the  main  features  of  this  approach  and  to  compare  the 
OVETV  algorithm  with  the  scalar  addition  and  scalar  multiplication  approaches.  In 
the  second  example,  identification  of  a  linear  time-varying  circuit  is  used  to  demon¬ 
strate  how  this  approach  can  be  applied  to  sampled-data  systems  with  quantization 
noise.  Finally,  we  show  how  this  approach  can  be  applied  to  actual  data  obtained 
from  a  crude  oil  distillation  column. 


3.2  First-order  Example 

Consider  the  following  time  varying  system. 


Uk 

=  -akVk-i  +  bkUk-i  +  Vk 

(3.1) 

ak 

—  cLk  -|-  0.0005 

(3.2) 

bk 

=  bk  +  omi 

(3.3) 
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where  Vk  is  a  uniformly  distributed  random  noise  between  [—0.1, 0.1].  For  this  simula¬ 
tion  Uk  is  a  uniformly  distributed  random  variable  between  [—1, 1].  Initial  conditions 
for  the  system  are  yo  =  0,  ao  =  0.5,  and  feo  =  l-O. 


To  apply  OVETV,  the  bounds  in  (2.3)  and  (2.4)  must  be  specified.  We  set  7^. 

iTlI 


-7  =0.1  and  Rk  =  CP  where  Cfc  = 


=  0.00118  and  I  is  the 


0.0005  0.001 

identity  matrix.  The  orientation  matrix  Pq  is  set  to  27,  and  the  center  6q  is  set  to 
the  origin. 


In  Figure  7  we  see  the  results  of  the  algorithm  at  one  time  step,  k  =  16.  We  begin 
with  the  ellipsoid  at  the  previous  time  step,  E15.  Applying  the  time  update  equations 
produces  the  ellipsoid  Exq.  The  measurements  U\q  and  yiQ  give  us  the  region  between 
two  hyperplanes,  Fiq.  Finally,  the  measurement  update  equations  are  used  to  find 

A 

the  minimum  volume  ellipsoid,  Eie,  which  bounds  the  intersection  of  Exe  and  Fxq. 
Note  that  the  -|-’s  represent  the  “true”  parameter  vector  at  ^  =  15  and  k  —  16. 

Results  of  the  simulation  for  A:  =  0  to  A:  =  300  are  shown  in  Figure  8  with  a*  at 
the  bottom  of  the  plot  and  hk  at  the  top.  The  upper  and  lower  bounds  are  found 
by  projecting  the  ellipsoids  onto  the  coordinate  axes  [29].  Simulation  results  for  the 
same  system  using  scalar  addition  and  scalar  multiplication  instead  of  the  optimal 
time  update  equations  are  shown  in  Figures  9  and  10,  respectively.  The  similarity 
of  these  plots  to  Figure  8  indicate  that,  for  this  system,  very  little  performance  is 
sacrificed  using  either  scalar  addition  or  scalar  multiplication. 

In  Figure  11,  the  volume  of  Ek  is  plotted  for  the  three  simulations.  Unlike  the 
linear  time  invariant  OVE  and  MOVE  algorithms,  OVETV  has  no  guarantee  of  mono- 


Parameters 


Parameters  Parameters 
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Figure  10:  Scalar  multiplication  results  for  First-Order  example 
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tonically  decreasing  volume.  For  this  example,  the  volume  of  Ek  using  scalar  addition 
is  almost  identical  to  the  volume  of  Ek  using  the  optimal  update  equations.  The  vol¬ 
ume  of  Ek  using  scalar  multiplication  is  larger,  but  the  difference  in  volume  is  small. 
To  ensure  that  the  volume  of  Ek  remains  “small”,  persistent  excitation  is  required. 
With  “less”  persistent  excitation,  the  difference  in  volume  between  the  algorithms 
would  probably  be  larger. 

The  normalized  center  estimate  error  for  all  three  simulations  is  plotted  in  Fig¬ 
ure  12.  This  error  is  calculated  according  to 

{0k-hfPk\^k-h) 

where  0k  and  Pk  parameterize  Ek,  and  0k  is  the  true  parameter  vector  [9].  If  the 
normalized  center  estimate  error  is  less  than  one,  then  the  true  parameter  vector,  0k, 
is  bounded  by  the  ellipsoid,  Ek',  this  consistency  is  guaranteed  by  all  three  algorithms. 


3.3  Linear  Time- Varying  Circuit 

Linear  circuits  with  time-varying  components  arise  in  various  applications.  For  ex¬ 
ample,  the  values  of  resistors,  capacitors,  and  inductors  all  vary  with  temperature. 
Resistances  and  capacitances  will  also  vary  with  adjustments  in  potentiometers  and 
variable  capacitors,  respectively. 

In  this  example,  we  consider  a  series  R-L-C  circuit  across  a  voltage  source  [30]. 
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Figure  12:  Normalized  center  estimate  error  for  First-Order  example 
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The  state  equations  for  this  system  are  given  as  follows: 


x{t) 

=  A{t)x(t) B{t)u{t), 

(3.4) 

2/(0 

=  D{t)x{t), 

(3.5) 

where 


1 


cW. 

m 


,  B(t)  = 

0 

1 

,  D(t)  = 

0  r{t) 

m  J 

(3.6) 


Xi{t)  is  the  voltage  across  the  capacitance  c{t),  X2{t)  is  the  current  through  the  in¬ 
ductance  /(t),  u{t)  is  the  voltage  source,  and  y{t)  is  the  voltage  measured  across  the 
resistance  r{t). 

To  evaluate  the  performance  of  our  algorithm,  it  is  necessary  that  we  discretize 
the  state  variable  equations  and  then  transform  them  to  an  ARX  structure.  Unfor¬ 
tunately,  for  linear  time-varying  systems,  in  general,  there  is  no  closed  form  solution 
for  discretizing  the  state  space  equations.  Consequently,  we  must  resort  to  numerical 
methods.  The  discretized  state  space  equations  have  the  form 


x{{k+l)T)  =  G{kT)x{kT)  +  H{kT)u{kT),  (3.7) 

y{kT)  =  D{kT)x{kT),  (3.8) 

where  T  is  the  sampling  time.  The  output  matrix  D{kT)  is  simply  the  output  matrix 
of  (3.5)  evaluated  at  time  kT.  To  find  the  state  matrix  G{kT)  and  the  input  matrix 
H{kT),  we  must  solve  numerically  for  the  state  transition  matrix,  using 


$(t,  kT)  =  A(t)#(t,  kT)  $(A;T,  kT)  =  1. 


(3.9) 
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The  matrices  G{kT)  and  H{kT)  can  then  be  evaluated  as  [31] 

G{kT)  =  $((A:  +  l)r,fcT),  (3.10) 

r{k+i)T  ^ 

H{kT)  =  mk  +  l)T,kT)j^^  ^-\T,kT)B{T)dT.  (3.11) 

Next,  we  need  to  transform  from  the  state  space  equations  to  an  ARX  structure. 
Because  the  system  matrices  in  (3.10)  and  (3.11)  are  time-varying,  it  is  not  possible 
to  find  the  pulse-transfer-function  in  the  usual  way  via  the  z  transform.  However, 
using  the  following  Theorem  (which  extends  a  concept  found  in  [30]  for  continuous 
time  systems)  we  are  able  to  transform  the  state  space  equations  of  (3.7)  and  (3.8) 
to  the  required  form. 


Theorem  5  Given  the  linear  time-varying  SISO  state  space  equation 

x{{k-\-l)T)  =  G{kT)x{kT)-i-H{kT)u{kT), 
y{kT)  =  D{kT)x{kT)-i-E{kT)u{kT), 


(3.12) 

(3.13) 


where  x{kT)  G  9?”  and  T  is  the  sampling  interval,  define  a  sequence  oflxn  matrices 
Lo{kT)  =  D{kT) 

L,(kT)  =  Li.i({k-H)T)G{kT)  i  = 

If  the  rank  of  the  matrix 


L{kT)  = 


Lo{kT) 

Li{kT) 


(S.15) 


L  L„-i(fcr)  J 

is  n  for  all  k,  then  y{kT)  satisfies  a  linear  nl^-order  difference  equation  of  the  form 


y{{k  -f-  n)r)  +  f)  ai{kT)y{{k  -F  n  -  i)T)  =  X)  fii{kT)u{{k  +  n-  i)T),  (3.16) 
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where 


an{kT)  an-i{kT)  •••  ai{kT) 


-Ln{kT)i-\kT), 


(3.17) 


'  MkT)  ^n-iikT)  •••  ^o{kT)]=MkT)-Ln{kT)l-\kT)W{kT),  (3.18) 


W{kT) 

Wn{kT) 


woo{kT)  0  •  •  •  0  O' 

wiQ{kT)  wii{kT)  : 

W{n-l)o{kT)  W^ri-l)l{kT)  •••  W(n-l){n-l){kT)  0_ 

Wno{kT)  WnlikT)  •••  Wn{n-l){kT)  Wnn{kT)  , 


(3.19) 

(3.20) 


and  the  Wij{kT)  are  given  recursively  by 

_  j  E{kT)  j  =  0 

-  I  X(,_i)o((fc  +  l)m(fcr)  i  =  l,2,---, 

Wij{kT)  =  u;(i_i)(j_i)((fc  +  l)r)  z  =  l,2, •••,  j  =  l,2, 

Proof:  The  proof  is  straight  forward.  From  equation  (3.13),  we  have 


(3.21) 

(3.22) 


y{kT)  =  D{kT)x{kT)  +  E{kT)u{kT) 

=  Lo{kT)x{kT)  +  woo{kT)u{kT),  .  (3.23) 

where  Lo{kT)  and  woo{kT)  are  as  defined  above.  Applying  the  shift  operator  q{-), 
where  q(x{kT))  =  x{{k  +  l)r),  to  equation  (3.23),  we  get 

y({k  +  1)T)  =  Lo{{k  +  l)T)x{{k  +  1)T)  +  woo{{k l)T)u{{k 1)T) 

=  Lo{{k  +  l)T){G{kT)ix{kT)  +  H{kT)u{kT)) 

+woo{{k  +  l)T)u{{k  +  l)r) 

=  Li{kT)x{kT)  +  wro{kT)u{kT)  +  wn{kT)u{{k  +  l)T).  (3.24) 
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Repeatedly  applying  the  shift  operator,  we  find  that  in  general, 

y((k  +  i)T)  =  Li(kT)x(kT)  +  ^  Wiju((k  +  j)T). 

j=o 

From  (3.25),  we  get  the  following  equations 


(3.25) 


y(kT)  =  L(kT)x(kT)-i-W(kT)u(kT) 
y((k-\-n)T)  =  Ln(kT)x(kT) -h  Wn(kT)u(kT). 


(3.26) 

(3.27) 


where  L(kT),  W(kT),  Ln(kT),  and  w„(kT)  are  as  defined  above,  and 


yikT)  = 

y{kT) 
y{{k  +  i)T) 

II 

Is" 

■  u{kT)  ■ 
«((A;-l-l)r) 

.  y{{k n  -  1)T)  . 

.  u{{k  -1-  n)T)  _ 

Solving  (3.26)  for  x(kT)  (L(kT)  is  invertible),  and  substituting  the  result  into  (3.27), 
we  get 


y((k  +  n)T)  -  Ln{kT)L  ^{kT)y{kT)  =  {wn{kT)  -  Ln{kT)L-'^{kT)W{kT))u{kT). 

(3.29) 

The  linear  n‘^-order  difference  equation  of  (3.16)  follows  directly  from  (3.29).  □ 
Note  that  the  condition  that  the  rank  of  L{kT)  be  n  guarantees  that  the  system 
IS  observable  (see  [30]  for  a  discussion  of  a  concept  in  continuous  time  referred  to 
as  instantaneously  ohservahle).  In  fact,  for  the  time-invariant  case  L  reduces  to  the 
observability  matrix,  i.e.,  U  =  [  J .  This  theorem  can 

be  extended  to  the  single-input,  multiple-output  (SIMO)  case  by  changing  L-\kT) 
in  the  definition  of  a,  and  A'  to  L+{kT),  where  (•)+  is  the  matrix  pseudoinverse  as 
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defined  in  [32].  In  the  next  chapter,  we  will  extend  this  theorem  to  multiple-input, 


single-output  (MISO)  systems. 

To  apply  Theorem  5  to  the  circuit’s  discretized  state  space  equations  m  (3.7)  and 


(3.8),  we  need  to  check  if 


Z(kT)  = 


0  diikT) 

d2((k  -F  l)T)g2i{kT)  d2{{k  -f  l)T)g22{kT),  J 


(3.30) 


is  invertible  where  G(-),  H{-),  and  D{-)  have  the  elements  gij{-),  and  d,(-), 

respectively.  Clearly,  L{kT)  will  be  invertible  provided  that  d2{kT)  and  g2i{kT)  are 
not  equal  to  zero  for  all  k.  This  will  be  the  case  for  the  system  we  consider. 

Applying  Theorem  5  to  (3.7)  and  (3.8),  we  get  the  following  difference  equation 

y{k  +  2)  +  ^  (3  3,1 

+^^^^^{g22{k)gn{k)  -  gi2{k)g2i{k))y{k)  =  (3-31) 

Mtt^lmib±ll(hi(k)gn{k)  -  h2{k)gn{k))u{k)  +  d2{k  -h  2)h2{k  -k  l)u{k  +  1) 

where  the  sampling  time,  T,  is  implied.  If  we  apply  the  inverse  shift  operator  q-'{-) 

twice  to  (3.31),  where  q-'{x{kT))  =  x{{h  -  l)r),  and  solve  for  y{k).  we  get  the 

following  equation 

=  P.32) 


where 


Ok  = 


<l>k 


■  -  2)52i(fc  -  1)  -  922{k  -  l)g2i{k  2)) 

~  ‘Z)g\\{k  -  2)  -  g\2{k  -  T)g2\{k  -  2)) 

—  2)g2i{k  —  2)  —  h2{k  —  2)gn{k  —  2)) 

g2l(k-2)  '■ 

iT 

[  — J/fc-l  —yk-2  Uk-1  Uk-2  ]  • 


(3.33) 


(3.34) 


The  discrete  time  input-output  equations  of  the  linear  time-varying  circuit  are 
given  by  equations  (3.9)-(3.11)  and  (3.32).  However,  because  the  OVETV  algorithm 
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is  implemented  on  a  digital  computer,  quantization  errors  in  the  D/A  and  A/D  con¬ 
verters  become  an  issue.  Consequently,  the  measured  input  Umk  and  the  measured 
output  Umk  differ  from  the  actual  input  into  the  system  Uk  and  the  actual  output  from 
the  system  yk-  They  are  related  by  the  following  equations 

U>k  '^mk  d”  '^qk  (3.35) 

Vk  ~  Vnik  l/qk  (3.36) 

where  Ugk  is  the  quantization  error  introduced  by  the  D/A  converter  and  ygk  is  the 
quantization  error  introduced  by  the  A/D  converter.  The  quantization  errors  are 
known  to  be  bounded  by 

(3.37) 

IVqkl  <  (3.38) 

where  the  quantization  level  Q  (for  the  A/D  and/or  D/A)  is  equal  to  FSAR 

is  the  full  scale  analog  range  of  the  converters,  and  nb  is  the  number  of  bits.  [31] 

In  general,  equation  (3.32)  can  be  written  as 

n  m 

yk  ^  ^  ^ikyk—i  “f"  ^  ^ik'^k—i^  (3.39) 

*=1  i=0 

After  substituting  equations  (3.35)  and  (3.36)  into  (3.39),  we  get 

n  m  n  m 

ymk  ^  V  ^ikym(k’-i)  H“  ^  ^  ^ik'^Tn(k’-i)  "f"  yqk  "t"  ^  y  ^ikyaik—t)  4”  ^  y  ^z7:^q(A:— z) *  (3.40) 

*=1  t=0  t=l  i=0 

If  we  let  Vk  =  yqk  +  ELi  ^ikyq{k-i)  +  ^ikUq^k-i)-,  equation  (3.40)  can  be  rewritten 
as 

n  m 

ymk  ^ikym{k—i)  “h  ^ik'^m(k—i)  4" 


(3.41) 
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which  is  precisely  the  ARX  structure  of  (2.1). 

To  apply  the  OVETV  algorithm,  we  need  a  bound  on  u*.  It  is  easy  to  show  that 

<  ^(1  + 1  l»«l)  +  %  E  (3-«) 

i=l  ^  i=0 

Consequently,  if  we  have  a  priori  bounds  on  the  parameter  values  (which  we  have 
already  assumed),  then  it  is  possible  to  derive  a  bound  on  Vk.  It  is  clear,  though,  that 
the  tighter  these  bounds  are,  the  better  the  OVETV  algorithm  will  perform. 

We  can  now  apply  OVETV  to  the  linear  time- varying  circuit.  We  will  consider 
the  case  where  r{t)  =  (1.0  -f-  0.5cos(f/10.0))0,  /(t)  =  l.OH,  and  c(t)  =  l.OF.  The 
sampling  time  T  is  selected  to  be  0.1  seconds.  Both  the  A/D  and  D/A  converter  are 
assumed  to  have  12  bit  resolution  and  a  FSAR  of  10  volts. 

The  OVETV  algorithm  is  initialized  as  follows.  Using  (3.42),  the  bounds  on  Vk 
are  set  at  7*  =  -7^.  =  0.005.  The  matrix,  Rfc,  which  bounds  the  parameter  variations 
is  set  at  Rk  =  CP  with  Cfc  =  0.001.  The  orientation  matrix,  Pq,  is  set  at  10/,  and  the 
ellipsoid  center,  6q,  is  set  at  the  origin.  For  this  simulation,  Umk^  the  input  supplied  by 
the  digital  computer,  is  selected  as  a  uniformly  distributed  random  variable  between 
[—5.0, 5.0]  volts. 

Results  of  the  simulation  are  shown  in  Figure  13  where  the  solid  lines  are  the 
true  parameters  calculated  from  (3.32)  and  the  dashed  lines  are  the  center  estimates 
obtained  from  the  OVETV  algorithm.  As  discussed  earlier,  with  this  approach  we 
are  more  interested  in  the  parameter  sets  than  the  center  estimates.  That  said, 
we  would  like  to  know  how  the  center  estimates  of  the  OVETV  algorithm  compare 
with  the  point  estimates  of  a  traditional  parameter  estimation  technique.  The  same 
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Table  2:  Performance  measures  for  OVETV-WRLS  comparison 


WRLS 

OVETV 

(Center  Estimate) 

£1 

£2 

0.0851 

0.0643 

0.0618 

0.0471 

simulation  was  run  using  weighted  recursive  least  squares  (WRLS).  The  WRLS  was 
tuned  empirically  to  minimize  the  following  performance  measures: 

1  ^  * 

=  '  =  (3-43) 

k=i 

where  N  is  the  number  of  samples,  Ok  is  the  parameter  estimate  at  time  k,  Oj,  is 
the  true  parameter  value  at  time  k,  and  ||  •  ||,-  is  a  vector  norm.  The  performance 
measures  resulting  from  a  simulation  with  a  forgetting  factor  of  0.73  and  an  initial 
covariance  matrix  of  10®  are  shown  in  Table  2.  The  performance  measures  for  the 
center  estimates  of  the  OVETV  algorithm  are  also  shown  in  Table  2.  The  center 
estimates  for  the  OVETV  algorithm  had  smaller  performance  measures  than  WRLS 
for  both  i  =  1  and  i  =  2. 

Despite  the  good  performance  of  the  center  estimates,  we  are  most  interested  in 
the  parameter  sets.  The  parameter  bound  for  bi  and  62  are  shown  in  Figures  14  and 
15,  respectively.  The  bounds  for  bi  are  especially  tight.  While  the  bounds  on  oj  and 
02  are  not  as  tight,  the  volumes  of  the  ellipsoids,  as  shown  in  Figure  16,  are  quite 


small. 


Parameters 


Ellipsoid  Volume 
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Figure  15:  Parameter  62  of  Linear  Circuit  example 
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Figure  16:  Ellipsoid  volume  for  Linear  Circuit  example 
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3.4  Crude  Oil  Distillation  Column 

As  a  final  example,  we  consider  a  crude  oil  distillation  column.  Our  identification 
exercise  for  this  example  will  be  conducted  on  actual  data  recorded  from  an  Amoco 
refinery  in  south  Texas.  The  goal  of  the  distillation  process  is  to  separate  a  mixture 
of  two  or  more  substances  into  its  various  components  with  a  desired  degree  of  purity. 
Because  different  pure  liquids  exhibit  different  volatilities,  the  application  and  removal 
of  heat  can  be  used  to  separate  the  components  [33]. 

While  the  distillation  process  is  highly  nonlinear,  within  a  given  operating  range 
it  can  effectively  be  modeled  as  a  linear  system.  However,  because  of  factors  such  as 
changes  in  raw  materials  or  production  levels  or  equipment  fouling,  the  parameters 
of  the  linear  model  tend  to  vary  with  time  [33].  Consequently,  for  control  purposes 
it  would  be  advantageous  to  track  these  parameters  as  they  change  with  time.  This 
information  could  be  then  be  utilized  to  retune  the  process  controllers  automatically, 
rather  than  manually,  as  is  often  done  in  practice. 

Several  data  sets  were  collected  under  open  loop  conditions  on  the  distillation 
column  shown  in  Figure  17.  The  results  from  one  test  are  shown  in  Figure  18.  The 
types  of  data  measured  are  given  in  Table  3  where  the  j/’s  are  outputs,  the  u’s  are 
inputs,  and  dl  is  a  disturbance.  This  is  clearly  a  multi-input,  multi-output  (MIMO) 
system. 

In  this  section,  we  will  look  at  the  SISO  relationship  between  the  heat  which 
is  being  applied  at  the  base  of  the  column,  u2,  and  the  top  section  temperature 
difference,  y2.  By  controlling  the  temperature  at  specific  locations  along  the  column. 
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Figure  17:  Distillation  Column 
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Figure  18:  Open  loop  data  for  Distillation  Column  example 


Table  3:  Distillation  column  data 


2/1 

overhead  drum  level 

y2 

top  section  temperature  difference  (degrees  F) 

2/3 

mid  section  temperature 

ul 

reflux 

u2 

heat  addition  to  base  of  tower 

dl 

feed  rate  to  tower 

3.5.  Summary 
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we  are  able  to  control  the  purity  of  the  product  compositions.  In  the  next  chapter, 
we  will  look  at  a  multi-input,  single-output  (MISO)  relationship  for  this  example. 

Based  on  knowledge  of  the  system  operation,  and  on  control  needs,  the  following 
model  structure  is  assumed 

Vk  =  -aij/fc-i  -  a2yk-2  +  huk-2  +  &3«fc-3,  (3.44) 

where  the  sampling  time  is  20  seconds.  The  OVETV  algorithm  is  initialized  as  follows: 
^k  =  “Xt  ~  ~  ~  ^-00001  >  Po  =  10/,  and  9o  =  0. 

Applying  the  OVETV  algorithm  to  this  data,  we  obtain  the  center  estimates  shown 
in  Figure  19.  While  we  cannot  compare  the  center  estimates  to  the  “true”  parameters 
(which  are  eflPectively  unknown),  we  can  observe  how  well  the  center  estimates  predict 
the  output  at  each  time  step.  In  Figure  20,  we  see  that  the  predicted  output  yk  =  Oj^k 
and  the  output  yk  are  almost  indistinguishable.  The  volumes  of  the  ellipsoids  are 
shown  in  Figure  21  and  the  parameter  bounds  at  298  minutes  are  shown  in  Figure  22. 

As  a  final  note,  for  the  bounds  specified  above  on  the  disturbance  vector,  Jk 
7^,  the  time-invariant  OVE  algorithm  was  inconsistent  and  failed. 


3.5  Summary 

In  this  chapter,  we  examined  three  examples.  The  first  example  was  a  first-order 
system  which  served  to  illustrate  the  key  features  of  the  OVETV  algorithm.  For  this 
example,  two  other  strategies,  scalar  addition  and  scalar  multiplication,  were  shown 
to  be  effective  alternatives  to  OVETV. 


Degrees  Fahrenheit  Parameters 
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Figure  19:  Center  estimates  for  Distillation  Column  example 


Figure  20:  Top  section  temperature  difference,  Distillation  Column  example 


imtifsi 


Figure  21:  Ellipsoid  volume  for  Distillation  Column  example 
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In  the  second  example,  we  explored  the  effects  of  sampling  and  quantization  on 
a  linear  time-varying  system,  a  linear  time-varying  circuit.  Theorem  5  was  given 
for  transforming  a  dicrete-time  linear  time-varying  state  space  equation  to  a  time- 
varying  difference  equation.  We  were  able  to  successfully  identify  the  system  using 
the  OVETV  algorithm.  In  fact,  while  the  goal  of  the  OVETV  algorithm  is  to  identify 
a  parameter  set,  the  center  estimate  of  the  OVETV  algorithm  was  able  to  track  the 
“true”  parameter  vector  better  than  the  “best”  estimate  using  WRLS. 

In  the  final  example,  we  applied  the  OVETV  algorithm  to  actual  data  obtained 
from  a  distillation  column.  Because  the  “true”  parameter  vector  was  effectively  un¬ 
known,  we  were  not  able  to  compare  it  with  the  parameter  sets  and  center  estimates. 
However,  the  predicted  output  based  on  the  center  estimate  matched  the  actual  out¬ 
put  very  closely. 


CHAPTER  IV 


Algorithm  Extensions 


4.1  Overview 

As  we  saw  in  the  last  chapter,  OVETV  is  a  very  effective  tool  for  PSE  of  time- varying 
systems.  However,  as  we  began  to  apply  OVETV,  our  experiences  challenged  us  to 
extend  and  improve  the  properties  of  this  algorithm.  In  this  chapter,  we  discuss 
several  of  these  modifications  and  adaptations  to  the  OVETV  algorithm. 

As  we  applied  OVETV  to  the  crude  oil  distillation  column  of  Section  3.4,  we 
discovered  that  each  of  the  outputs  was  actually  affected  by  more  than  one  input. 
Clearly,  it  would  be  good  to  identify  this  relationship  directly  in  the  OVETV  frame¬ 
work.  In  Section  4.2,  we  discuss  how  the  OVE  and  OVETV  algorithms  can  be  used 
to  identify  multiple-input,  single-output  (MISO)  systems.  We  begin  by  showing  how 
MISO  linear  time-invariant  systems  can  be  placed  within  the  OVE  framework.  We 
also  extend  Theorem  5  of  Section  3.3  to  show  how  certain  MISO  linear  time-varying 
systems  can  be  placed  within  the  MISO  framework.  Finally,  we  use  the  results  of  this 
section  to  identify  a  MISO  relationship  for  the  distillation  column  example. 

Often,  as  was  the  case  in  the  linear  time-varying  circuit  example  of  Section  3.3, 
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the  time- varying  parameters  that  we  axe  trying  to  identify  are  dependent  on  a  smaller 
set  of  physical  time-varying  parameters.  If  we  can  exploit  this  dependence,  we  may 
be  able  to  simplify  the  problem  and  reduce  the  number  of  computations  required  to 
implement  the  algorithm.  Adaptations  of  the  OVETV  algorithm  are  developed  in 
Section  4.3  to  utilize  knowledge  of  these  dependencies. 

When  applying  the  OVETV  algorithm  in  practice,  the  orientation  matrix,  Pk,  can 
become  non-positive  definite.  Theoretically  this  is  impossible,  geometrically  it  does 
not  make  since,  and  practically  it  causes  the  algorithm  to  fail.  Numerical  round¬ 
off  errors  are  causing  the  algorithm  to  bomb.  To  solve  this  problem,  we  apply  an 
approach  that  has  been  used  successfully  for  stochastic  estimation.  A  “square-root” 
algorithm  is  developed  for  implementing  OVE.  Furthermore,  it  is  shown  how  the 
square-root  approach  can  be  extended  to  time- varying  systems. 

The  type  of  input  into  the  system  has  a  very  significant  effect  on  the  “quality” 
of  estimates  produced  from  an  identification  scheme.  A  “calculated”  synthesis  of 
the  system  input  is  explored  in  Section  4.5.  In  particular,  we  investigate  an  OVE- 
based  input  synthesis  procedure  (OVE-ISP)  which  was  developed  in  [29]  and  [21]. 
We  show  how  the  scope  of  OVE-ISP  can  be  extended  to  handle  systems  with  known 
transportation  lag  such  at  the  distillation  column  in  Section  3.4.  We  also  demonstrate 
its  application  to  time-varying  systems  and  an  important  “stabilizing”  property  that 
other  synthesis  strategies  do  not  appear  to  possess. 


4.2.  Multiple-Input,  Single-Output  Systems 
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4.2  Multiple- Input,  Single-Output  Systems 

In  many  real  systems,  the  output  of  the  system  is  effected  by  more  than  one  input. 
Consequently,  we  would  like  to  extend  the  OVE  and  OVETV  algorithms  to  handle 
MISO  systems.  In  this  section  we  will  show  how  this  can  be  done  for  linear  time- 
invariant  systems.  Then  we  will  show  how  this  can  be  done  for  a  class  of  linear  time- 
varying  systems.  Finally,  we  will  apply  the  MISO  framework  to  the  crude  distillation 
column  example  which  we  explored  in  Section  3.4. 


4.2.1  Linear  Time-Invariant  Case 


Consider  the  linear  time-invariant  MISO  state  space  equation 


x{{k-\-l)T)  =  Gx{kT)  +  f2HnU,{kT),  (4.1) 

«=1 

y{kT)  =  Dx{kT)  +  j2E.u,,{kT),  (4.2) 

K=1 


where  x{kT)  €  u^ikT)  6  {G,  are  matrices  of  appropriate  dimen¬ 

sions,  and  T  is  the  sampling  interval.  Taking  the  Z  transform  of  (4.1)  and  (4.2) 
gives 


zX{z)  =  GX{z)  +  f2HMz),  (4.3) 

«=1 

Y{z)  =  DX{z)  +  f2EnU.{^)- 

K=1 


(4.4) 


74 


CHAPTER  IV.  Algorithm  Extensions 


Solving  (4.3)  for  X{z),  substituting  it  into  (4.4),  using  a  classical  representation  for 
the  matrix  inverse,  and  rearranging  terms  results  in 

det{zl  -  G)Y{z)  =  f^{i)adj(z/  -  G)H,  +  det(z7  -  G)E,]U,{z),  (4.5) 

K=1 

where  det(-)  and  adj(-)  return  the  determinant  and  adjoint  of  a  matrix,  respectively. 
If  we  multiply  (4.5)  by  take  the  inverse  Z  transform,  and  solve  for  y{kT),  it  is 
easy  to  see  that  we  get  an  equation  of  the  form 

vikT)  =  -  f:  aiviik  -  i)T)  +  ^212  -  i)n  (4.6) 

2  =  1  «=1  2=0 

Equation  (4.6)  can  be  rewritten  as 

yk  =  9'^<f>k  (4.7) 

where 

r  1  ^ 

e  =  [  ai  •  •  •  ct„  boi  ■■■  bni  .  6om  ' '  •  J  ,  (4.8) 

=  [—y{k  —  l)  ■■■  —y{k  —  n)  Uo{k)  ■■■  uo{k  —  n) 

.  u^i{k)  •••  u^{k-n)Y-  (4.9) 

Assuming  that  the  noise  or  uncertainty  entering  the  system  as  {yk  =  9"^ (j>k  +  Vk)  is 
bounded,  the  OVE  algorithm  can  be  directly  applied  to  the  system  described  by  (4.7). 

4.2.2  Linear  Time- Varying  Case 

The  goal  of  the  dissertation  has  been  to  identify  the  parameters  of  time-varying 
systems.  Clearly,  we  would  like  to  be  able  to  do  this  for  MISO  systems  as  well.  To 
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see  how  MISO  systems  might  fall  into  our  time- varying  ARX  framework,  we  will 
extend  Theorem  5  of  Section  3.3  to  MISO  systems. 


Theorem  6  Given  the  linear  time-varying  MISO  state  space  equation 


x{{k  +  l)T)  =  G{kT)x{kT)-lJ2H^(^T)u,,i^T), 


y{kT)  =  D{kT)x{kT)-\-j^E,,{kT)u,,{kT), 


(4.10) 


(4.11) 


where  x{kT)  G  u,,{kT)  €  and  T  is  the  sampling  interval,  define  a  sequence  of 


1  X  n  matrices 


Lo{kT)  =  D{kT) 

Li{kT)  =  Li^y{{k  +  l)T)G{kT)  i  = 


(4.12) 


If  the  rank  of  the  matrix 


Lo{kT) 

HkT) 


L{kT)=  .  (4-13) 

.  T„_i(fcT)  . 

s  n  for  all  k,  then  y{kT)  satisfies  a  linear  n^^-order  difference  equation  of  the  form 

y{{k  -t-  n)T)  -\-f^ai{kT)y{{k  +  n  -  i)T)  =  I^UkT)un{{k  +  n  -  0^),  (4.14) 


KZ=1  i—0 


where 


(4.  IS) 


'  a„{kT)  an-i(kT)  a,(tT)  ]  = -l„(l:T)l-\l:T),  (4.i 

\MkT)  ■■■  Mf‘T)\=<l>„(kT)-L.,(hT)L-'{kT)WJkT), 


(4.16) 
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W,{kT) 

Wnn{kT) 


WoonikT)  0  •  •  •  0  O' 

Wio^{kT)  Wii^{kT)  '  • .  : 

.  '^(n— 1)Ok(^2'')  ^(n— l)lfv(^2^)  ■  ■  '  ^(n— l)(n— 0 


WnOnikT)  Wnln{kT) 


Wn{n-l)ti{kT)  W 

nnn  m  , 


(118) 


and  the  WijK{kT)  are  given  recursively  by 


wojf,{kT)  = 


E.{kT) 


i  =  o 


L(y_a)o((fc  +  l)m«(fcr)  i  =  l,2, 


w 


(kT)  =  Ki(, + l)r)  j  =  1,2, 


(4.19) 

(4.SO) 


and  K  =  1,2,  -  •  • ,  g. 

Proof:  The  proof  is  a  direct  extension  of  the  proof  of  Theorem  5  in  Section  3.3  and 
is  therefore  omitted.  Below  we  will  show  how  this  theorem  can  be  used  to  transform 
a  MISO  system,  which  can  be  realized  via  linear  time-varying  state  space  equations, 
to  a  framework  for  which  the  OVETV  algorithm  is  applicable.  □ 

If  we  apply  the  inverse  shift  operator,  to  (4.14)  n  times  and  solve  for  y{kT), 

we  get  the  following 

y{kT)  =  -f2  ai{kT)y{ik  -  i)T)  +  ^211  bi.{kT)u,{{k  -  i)T),  (4.21) 

«=1 

where  ai{kT)  =  ai{{k  —  n)T)  and  bi^ikT)  =  /di^dk  —  n)T).  Equation  (4.21)  can  be 
rewritten  as 

yk  =  eUk  (4.22) 

where 


ai{kT)  ■■■  an{kT)  boi{kT)  •••  6„i(fcT) 
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bou.{kT)  ■  ■  ■  bn^{kT) 


(t>k  = 


—y{k  —  l)  •••  —y[k  —  n)  Uo{k)  •••  uo{k  —  n) 
.  u^{k)  ■  •  •  Ufj,{k  -  n)  ]  . 


(4.23) 


(4.24) 


If  the  noise  or  uncertainty  entering  the  system  as  (?/*  —  Sk^k  +  Vk)  is  bounded,  and  if 
we  can  find  a  bound  on  how  much  the  parameters  change  during  any  time  step,  then 
the  OVETV  algorithm  can  be  directly  applied  to  the  system  described  by  (4.22). 


Remark  1  Before  giving  a  MISO  example,  this  seems  to  be  an  appropriate  place 
to  discuss  the  more  general  multiple-input,  multiple-output  (MIMO)  case.  First, 
Theorem  6  can  be  extended  to  the  MIMO  case  by  changing  L~^{kT)  in  the  definition 
of  ai  and  to  L'^{kT),  where  (•)■•■  is  the  matrix  pseudoinverse  as  defined  in  [32]. 
The  system  can  then  be  transformed  into  the  full  polynomial  form  as  defined  in  [34]. 
In  the  MIMO  case,  the  system  disturbance,  Vk,  is  a  vector.  If  bounds  are  known  on 
each  of  the  components  of  Vk,  then  the  problem  can  be  composed  into  a  separate 
problem  for  each  output  for  which  the  MISO  results  discussed  above  directly  apply 
[5].  When  a  bound  exists  on  the  norm  of  Ufc,  rather  than  on  each  of  the  components, 
the  problem  becomes  more  difficult.  For  extensions  of  the  OBE  algorithm  to  the  case 
where  a  norm  bounds  exists  on  Vk,  see  [5]. 
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4.2.3  Crude  Oil  Distillation  Column  Example 

In  Section  3.4,  we  estimated  the  parameters  of  the  difference  equation  relating  the 
heat  which  is  being  applied  at  the  base  of  the  distillation  column,  u2,  to  the  top  section 
temperature  difference,  y2.  In  actuality,  the  top  section  temperature  difference  is  also 
affected  by  the  reflux,  ul,  and  the  (disturbance)  feed  rate  to  the  tower,  dl.  Clearly, 
we  are  interested  in  the  MISO  relationship  between  the  output,  y2,  and  the  inputs, 
ul,  u2,  and  dl. 

Based  on  knowledge  of  the  system  operation,  and  on  control  needs,  the  following 
model  structure  is  assumed: 

y2{kT)  =  -j2<kT)y2{{k-i)T)  +  j2bn(kT)uliik-i)T) 

i=l  i=2 

5  5 

+  E  bi2{kT)u2{{k  -  i)T)  +  E  bi^{kT)dl{{k  -  i)T)  (4.25) 

1=2  4=2 

where  the  sampling  time  is  20  seconds.  The  OVETV  algorithm  is  initialized  as  follows: 
Ta:  =  “Tfc  =  0.15,  Rk  =  (kl  with  Cfc  =  0.0000006,  Pq  =  107,  and  =  0. 

Applying  the  OVETV  algorithm  to  the  data  shown  in  Section  3.4,  we  get  the 
results  shown  in  Figure  23  for  two  parameters,  a\  and  622-  The  center  estimates  and 
parameter  bounds  for  all  the  parameters  at  298  minutes  are  shown  in  Figure  24.  As 
in  the  SISO  case,  we  cannot  compare  the  estimated  parameter  sets  with  the  “true” 
parameters  because  they  are  effectively  unknown.  Also,  as  in  the  SISO  case,  we 
can  observe  how  well  the  center  estimates  predict  the  output  at  each  time  step.  As 
would  be  expected,  the  predicted  output,  yk  =  Ok^ki  when  we  consider  multiple 
inputs,  matches  the  actual  output,  yk,  better  than  the  predicted  output  when  we 
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Minutes 


Figure  23:  Parameters  gj  and  622  of  MISO  Distillation  Column  example 

only  consider  a  single  input.  In  fact,  the  absolute  error,  \yk  —  yk\i  is  bounded  by 
7^  =  —7^  =  0.15  for  the  MISO  case.  The  volumes  of  the  ellipsoids  are  shown  in 
Figure  25. 

4.3  Dependent  Parameter  Variations 

In  this  section,  we  will  present  an  algorithm  which  exploits  information  regarding 
dependency  in  parameter  variations  to  reduce  the  computational  complexity  of  the 
OVETV  algorithm.  In  the  problem  statement  given  in  Chapter  II,  it  was  assumed 
that  jRfc,  the  matrix  bound  on  the  change  in  the  parameter  vector  during  each  time 


4.3.  Dependent  Parameter  Variations 


81 


step,  was  positive  definite.  While  the  theory  of  [15]  still  holds  if  Rk  G  is  positive 
semi- definite,  the  bound  on  Wk  must  be  defined  in  terms  of  its  support  function  [26] 
instead  of  (2.4),  that  is, 

{wk  e  di’’  :  rjjwk  <  ^/rijR^;  rjr  G  (4.26) 

The  case  where  Rk  is  positive  semi-definite  is  important  because  often  there  are  time- 
invariant  parameters,  or  the  variation  in  one  parameter  is  dependent  on  variations  in 
other  parameters. 

Furthermore,  if  it  is  known  that  the  rank  of  Rk  is  always  less  than  r,  we  can  use 
this  information  to  reduce  the  number  of  computations  required  for  the  time  update 
equations.  If  the  rank  of  Rk  is  always  s  <  r,  then  (2.2)  can  be  rewritten  as 

^fc+i  =  +  BkWki  (4-27) 

where  Wk  €  3^®  is  assumed  to  satisfy 

wlRk^Wk  <  1,  (4.28) 

and  the  matrix,  Rk  G  is  symmetric,  positive  definite  and  known  at  each  time  k. 
The  bound  in  (4.28)  can  also  be  rewritten  in  terms  of  its  support  function  as 

{wk  €  3?*  :  Tjjwk  <  ^/r)jRkris-,ris  €  3ft®}.  (4.29) 

Using  Wk  =  BkWk,  the  bound  in  (4.29),  and  the  definition  of  a  support  function,  the 
bound  in  (4.26)  is  given  by 


{wk  e  :  rjjwk  <  \J rij BkRkBlr}r\ rjr  G  3fJ’'}. 


(4.30) 
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Consequently,  the  generalized  eigenvalue  problem  in  (2.79)  can  be  rewritten 

Bk-\Rk-iBl_.^Xj  =  XjPk-iXj.  (4.31) 

Solving  (4.31)  for  each  Aj,  j  G  [l,r]  is  equivalent  to  solving 

det(APfc_i  -  Bk-iRk-xBl_^)  =  0  (4.32) 

for  A.  Using  properties  of  the  determinant  [35],  it  is  easy  to  show  that 


^^i{XPk-i  - Bk-iRk-iBl_Pj  =  det(.Rfc_i)det(A-i)A'-Met(A.R^.^i 

(4.33) 

Using  (4.31)  and  (4.33),  the  revised  time  update  equations  are  given  in  the  fol¬ 
lowing  procedure: 

1.  Solve  the  generalized  eigenvalue  problem 

^k—l^k—l^k-l^j  —  \jRf._yXj.  (4.34) 

for  each  Ay,  j  e  [l,^]. 


2.  Solve 


for  the  unique  p  >  0. 


E 


j=i 


1 

Ay  -fp 


{s  —  r)p  +  s 
p(p  -I- 1) 


(4.35) 


3.  Then,  the  solution  to  (2.11)  is  given  by 


=  h-\  (4.36) 

A  =  (p  +  l)Pi-i  +  (p“^  +  l)BkRk-iBl  (4.37) 
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Note  that  the  generalized  eigenvalue  problem  in  step  1  has  been  reduced  from 
rth  order  to  sth  order.  Furthermore,  the  (r  +  1)  order  polynomial  in  step  2  has 
been  replaced  with  a  (s  +  1)  order  polynomial.  If  s  <  r,  the  reduction  in  number 
of  computations  will  be  substantial.  Note  that  the  inverse  of  Rk-i  can  often  be 
calculated  a  priori.  For  the  particular  case  when  s  =  1,  see  [15]. 


4.4  Square-Root  Implementation 

One  problem  when  implementing  many  recursive  estimation  algorithms  is  that  the 
covariance  matrix  in  stochastic  algorithms,  such  as  the  Kalman  filter  and  recursive 
least  squares,  and  the  ellipsoid  matrix  in  bounded  input  algorithms,  such  as  OVE  and 
OVETV,  tend  to  loose  their  positive  definite  (PD)  property  due  to  finite  precision 
calculations.  This  may  cause  the  algorithms  to  “fail.”  One  solution  to  this  problem 
is  to  factor  the  PD  matrix  as  P  =  SS'^  and  then  update  S  instead  of  P;  this  would 
guarantee  the  PD  property  of  P.  There  have  been  many  square-root  methods  de¬ 
veloped  for  updating  the  covariance  matrix  of  stochastic  algorithms.  For  example, 
see  [36],  [37],  [38].  While  some  of  the  algorithms  are  easily  amenable  to  the  OVE 
algorithm,  others  are  not. 

In  this  section,  we  will  show  how  Potter’s  implementation  of  the  Kalman  filter  [36] 
can  be  modified  to  work  with  the  OVE  algorithm.  We  will  also  compare  the  number  of 
floating  point  operations  using  the  direct  implementation  of  the  OVE  algorithm  with 
those  using  the  new  square-root  algorithm.  Finally,  we  will  show  how  the  square-root 
algorithm  can  be  extended  to  the  time-varying  case. 
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4.4.1  Linear  Time-Invariant  Case 


Recall  from  Section  2.3  that  the  measurement  update  equations  for  the  center  esti¬ 
mate,  Ok,  and  the  ellipsoid  matrix,  Pk,  are  given  by 


TkPk<t>k 


Ok  —  Ok  -\- 

Pk  =  4  A  +  (cTfc  —  h) 


Pk<l>k<l>kPk 

<i>lPk<t>k  ’ 


(4.38) 

(4.39) 


where  Tk,  8k,  and  Ok  are  scalars  defined  in  Section  2.3  and  <j)k  is  the  regression  vector. 
For  the  time-invariant  case  Ok  =  Ok-i  and  Pk  =  Pk-i- 

Following  an  outline  similar  to  the  development  of  Potter’s  algorithm  for  the 
Kalman  filter  in  [36],  we  will  develop  a  square- root  algorithm  for  the  OVE  algorithm 
We  begin  by  setting 


Pk  =  SkSl 

Pk  =  SkSl 

Substituting  (4.40)  and  (4.41)  into  (4.39),  we  get 

Pk  =  SkSk  =  SkSkSk  +  (cTfc  -  8k)SkSl^kh^'kSkSl 

=  Sk{8kl  +  {cTk  -  8k)\k'cokXZ!k)Sk , 

where  Wk  =  Sld'k  and  At  =  —7! — . 

From  (4.42)  it  is  easy  to  see  that  we  can  define 


(4.40) 

(4.41) 


(4.42) 


Sk  =  Sk{8kl  +  {ok  -  6k)\k^k'^kY^^ 
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=  8l''^Sk{I  -  (1  -  (4.43) 

where  we  know  that  6k  is  positive.  Next,  we  use  the  fact  that  (7^(1  —  ^)XkWkwl) 
can  be  factored  as 

/  -  (1  -  ^)AfctJ7fcWfc  =  (/  -  (4.44) 

Ok 

where  flk  is  the  root  of 

zu^Wkl^l -2j3k  +  {1  - —0.  (4.45) 

The  roots  of  (4.45)  are  given  by  /?*  =  (1  ±  Using  the  root  0k  =  {1  +  \/^)Afc 

and  the  factorization  in  (4.44),  (4.43)  can  be  written  as 

Si  =  + 

=  -  (1  + 

where  Lk  =  Sk'^^k- 

Equation  (4.46)  is  key  step  in  the  square-root  algorithm  for  the  OVE  measurement 
update.  The  complete  algorithm  (excluding  the  calculation  of  the  scalars  r^,  4  and 
cTfe  which  can  be  found  in  Section  2.3)  is  shown  in  Table  4.  Also  shown  in  Table  4 
is  a  tabulation  of  the  number  of  floating  point  operation  required  for  each  step. 
For  comparison  purposes,  Table  5  details  the  steps  and  corresponding  floating  point 
operations  for  a  direct  implementation  of  the  OVE  measurement  update.  In  our 
floating  point  counts  for  the  direct  implementation,  we  take  advantage  of  the  fact  that 
we  only  have  to  calculate  the  upper  triangular  part  of  Pk  because  it  is  symmetric. 
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Table  4:  Square-root  algorithm  and  corresponding  floating  point  operations 


Steps 

Floating  Point  Operations 

+ 

X 

V" 

Wk  =  Skh 

'n? 

Xk  —  — T - 

n 

n 

\T  = 

1 

1 

Lk  =  Sk'^k 

n? 

n? 

+  LkSk 

n 

n 

sT  =  Vh 

1 

1 

1 

1 

Sk  =  ^y'^Sk  -  {ukLk)wl 

n? 

2n^  +  n 

Total 

372^  -|-  272  -|-  1 

4n^  -}-  3n  -f  2 

3 

1 

As  can  be  seen  in  Table  4  and  Table  5,  the  number  of  floating  point  operations  is 
of  the  same  order  for  each  algorithm.  The  number  of  floating  point  operations  using 
the  square-root  algorithm  will  be  more  than  the  number  of  floating  point  operations 
using  the  direct  implementation.  In  most  cases,  however,  this  increase  in  number  of 
computations  will  be  less  important  than  the  increase  in  numerical  stability  which  is 
gained  by  using  the  square-root  algorithm. 

Finally,  we  need  to  be  concerned  with  how  the  algorithm  is  initialized.  That  is, 
how  does  one  find  Sq  assuming  that  Pq  is  specified.  Often,  when  little  is  known  about 
the  plant,  Pq  is  initialized  to  ^7  where  I  is  the  identity  matrix  and  ^  is  a  very  large 
positive  scalar.  In  this  case.  So  can  simply  be  set  to  a/^/.  If,  on  the  other  hand,  more 
is  known  about  the  system  and  Pq  has  a  more  complicated  structure,  then  Pq  can 
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Table  5:  Direct  implementation  and  corresponding  floating  point  operations 


Steps 

Floating  Point  Operations 

+ 

X 

— 

~  Pk^k 

in? 

n 

n 

1 

1 

St  =  4 

n 

n  -h  1 

L]t  —  ^k')^k 

1 

n-\-2 

Pk  =  4  A  +  Lkwl 

in2  +  in 

-t-  n 

Total 

2n^  -|-  4n  -|-  3 

1 

1 

easily  be  factored  as  SqSq 


using  the  Cholesky  decomposition  [32]. 


4.4.2  Linear  Time- Varying  Case 

Because  of  the  improved  numerical  stability  of  square-root  algorithms,  we  would  like 
to  extend  this  concept  to  the  time-varying  case.  Unfortunately,  the  optimal  time 
update  equations  of  the  OVETV  algorithm  are  not  amenable  to  being  rewritten  in 
square-root  form.  However,  the  two  alternative  algorithms  discussed  in  Section  2.6, 
scalar  addition  and  scalar  multiplication,  are  very  amenable  to  being  rewritten  in 
square-root  form. 
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Scalar  Addition 


Recall  from  Section  2.6  that  the  basic  idea  of  scalar  addition  was  to  reduce  the  number 
of  calculations  by  constraining  the  time  update  ellipsoid,  Ek,  to  be  an  ellipsoid  with 
the  same  center  and  orientation  as  Ek-iy  and  with  semi-axis  lengths  which  are  found 
by  adding  a  scalar,  rj,  to  the  semi-axis  lengths  of  Ek-i.  Rk  is  also  constrained  to  be 


ai- 

From  Theorem  3,  we  see  that  this  is  equivalent  to  constraining  Pk  to  be  of  the 
form 


(4.47) 


where  a  SVD  of  Pk-\  is  given  by 


A-i  =  Vk-,Ak-,Vl,, 


(4.48) 


14-1  is  orthogonal,  and  Afe_i  =  diag(Ai . . .  A^). 

To  show  how  scalar  addition  can  be  cast  in  a  square-root  framework  we  take  a 
SVDof  5fc_a, 

Sk-\  —  (4.49) 

where  14_i  and  Uk-i  are  orthogonal  and  Ek-i  —  diag(cri . . .  (7r)-  Substituting  (4.49) 
into  (4.40),  we  get 


Pk-i  =  Sk-iSl_, 

=  Vk.^Ek-iUl,Uk-iEk-iVl, 

=  14-iStilf_i. 


(4.50) 
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Comparing  (4.48)  and  (4.50),  we  see  that  Vk-i  and  Afc_i  can  be  chosen  such  that 
Vk-i  =  Vfc-i  and  =  Sfc_i. 

Let 

5fc  =  T4_i(S,_i+777)Ci-  (4-51) 

Substituting  (4.51)  into  (4.41),  we  get 


ft  =  hsj 

=  +  r,l)Ut,Uk-i(T.k-i  + 

=  (4-52) 


1  /2 

Comparing  (4.47)  and  (4.52),  and  using  the  fact  that  14_i  =  14_i  and  A^_j  =  S^-i, 
we  see  that  (4.51)  provides  us  with  the  square-root  update  for  scalar  addition  where 
the  SVD  of  Sk-i  is  given  by  (4.49). 

However,  we  still  must  show  that  the  value  of  rj  which  minimizes  the  volume  of 
Ek  such  that  Ek  3  Gk  where  Gk  is  given  by  (2.26)  can  be  calculated  directly  from 
Sk-i-  From  Theorem  3,  we  know  that  optimum  value  for  t]  is  given  by 


V  =P6+P7- 


2aa 

3(u  “|-  o) 


(4.53) 


where 


a=max\/Xj,  a=  min  JXj,  (4.54) 

je[i,r]V  je[i,r]V 

and  pq  and  pr  are  defined  in  Theorem  3  and  are  functions  of  a,  a,  and  (k-\  ■  Therefore, 
if  we  can  define  a  and  a  in  terms  of  Sk-i  we  have  achieved  our  goal.  In  fact, 

a  =  max  <t,-,  a  =  min  cr,-,  (4.55) 
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where  the  aj  are  the  singular  values  given  by  (4.49). 

Scalar  Multiplication 

In  the  scalar  multiplication  strategy,  Ek  is  constrained  to  be  an  ellipsoid  with  the 
same  center  and  orientation  as  Ek-i',  the  semi-axis  lengths  are  found  by  multiplying 
the  semi-axes  of  Ek-i  with  a  scalar,  r].  From  Theorem  4,  we  see  that  this  is  equivalent 
to  constraining  Pk  to  be  of  the  form 

Pk  =  rj^Pk-i.  (4.56) 

It  is  easy  to  see  that  this  is  equivalent  to 

Sk  =  r]Sk-i,  (4.57) 

in  our  square-root  framework. 

We  still  must  show  that  the  value  of  rj,  which  minimizes  the  volume  of  Ek  such 
that  Ek  D  Gk  where  Gk  is  given  by  (2.10)  {Rk  is  not  constrained  to  Ckl),  can  be 
calculated  directly  from  Sk-i-  From  Theorem  4,  we  know  that  optimum  value  for  rj 
is  given  by 

ri  =  lA\/i  (4.58) 

where  A  is  the  maximum  generalized  eigenvalue,  Aj,  which  satisfies 


Rk — 1  ^  j  —  ^jPk — 1  ^  j  * 


(4.59) 
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Since  Rk-i  is  symmetric  and  positive  semi-definite,  it  can  be  factored  as  Tk-iT^_^. 
Therefore,  (4.59)  can  be  written  as 

Tk-\Tk-\Xj  =  XjSk-iS^-iXj.  (4.60) 

This  problem  can  be  solved  directly  from  Tk-i  and  Sk-i  using  the  generalized 
singular  value  decomposition  [32].  Given  Tk-i  and  Sk-i,  there  exist  orthogonal  T4-i 
and  Uk-i  and  invertible  Xk-i  such  that 

Vl,Tk-iXk-i  =  Ck-k  (4.61) 

Ul,Sk-iXk-i  =  Dk-i,  (4.62) 

where  Ck-i  =  diag(ci . . .  c^)  and  Dk-i  -  diag(di . . .  dr).  The  generalized  singular 
values  are  given  by  crj  =  cj/sj  where  j  =  The  generalized  singular  values,  aj, 

are  also  equal  to  the  square- root  of  the  generalized  eigenvalues,  Xj.  Therefore,  the 
optimum  value  of  rj  is  given  by 

T}  =  1  a  (4.63) 

where  a  is  the  maximum  generalized  singular  value.  For  a  discussion  of  how  the 

generalized  singular  values  can  be  calculated,  see  [32]. 

Finally,  some  discussion  of  how  Rk  can  be  factored  as  is  in  order.  Clearly 

if  Rk  =  Ckl,  then  Tk  =  (kl-  Actually,  if  Rk  =  (ll-,  then  this  problem  can  be  solved 
using  the  standard  singular  value  decomposition.  If  on  the  other  hand,  Rk  has  a  more 
complicated  structure,  then  Rk  can  be  factored  using  the  Cholesky  decomposition  [32]. 
If  Rk  does  not  vary  with  time,  then  this  can  be  done  once  before  the  identification 
process  begins. 
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4.5  Input  Synthesis 

The  input  into  a  system  during  an  identification  experiment  greatly  affects  the  “qual¬ 
ity”  of  the  estimates  produced  by  an  identification  algorithm.  Often  the  experiment 
design  has  some  freedom  on  how  this  input  is  chosen.  We  will  call  this  design  of 
the  system  input,  which  is  intended  to  improve  the  identification  experiment,  “input 
synthesis.” 

In  this  section,  we  will  consider  an  input  synthesis  strategy  developed  to  work 
with  the  OVE  algorithm  [11].  The  goal  of  this  OVE-based  input  synthesis  procedure 
is  to  drive  the  system  such  that  the  regression  vector,  is  aligned  parallel  to  the 
ellipsoid  axis  of  greatest  length.  This  direction  is  chosen  because  it  was  found  in  [29] 
that,  when  using  the  OVE  algorithm,  the  greatest  reduction  in  volume  occurs  along 
the  direction  which  is  parallel  to 

The  OVE  Input  Synthesis  Procedure  (OVE-ISP)  was  originally  developed  in  [29] 
for  systems  with  a  single  delay,  /  =  1.  The  algorithm  is  based  on  a  system  model 
which  is  realized  from  the  ellipsoid  center  estimate  at  time  ko,  Oko^  for  which  <f)k  is 
the  state  vector.  In  [39],  DeVilbiss  and  Yurkovich  modified  the  algorithm  to  handle 
systems  where  the  delay  index,  I,  was  either  0  or  1.  For  improved  numerical  robust¬ 
ness,  they  used  a  balanced  minimal  realization  instead  of  the  canonical  realization 
discussed  above.  However,  by  using  a  minimal  realization,  the  new  strategy  required 
an  additional  “observation”  stage  which  was  not  needed  with  the  canonical  realiza¬ 
tion.  In  [11],  DeVilbiss  and  Yurkovich  returned  to  the  canonical  realization  because 
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it  was  found  that  the  input  signal  level  during  the  “observation”  stage  was  large  in 
magnitude  compared  to  the  final  stage  of  the  input  synthesis  strategy. 

This  input  synthesis  strategy  is  investigated  in  this  dissertation  for  several  rea¬ 
sons.  First,  the  authors  of  [11]  and  the  author  of  this  dissertation  feel  that  one  of 
the  most  significant  applications  of  OVE-ISP  may  be  for  robust  adaptive  control  of 
time-varying  systems.  The  reasoning  is  as  follows.  As  has  been  said,  when  applying 
PSE  to  robust  control  or  robust  adaptive  control,  there  is  a  trade  off  between  set 
size  and  system  performance.  In  [39],  the  OVE-ISP  algorithm  demonstrated  signifi¬ 
cantly  greater  ellipsoid  volume  reduction  than  a  pseudorandom  white  noise  sequence 
of  equal  expected  energy  during  the  “transient  phase”  of  the  identification  exper¬ 
iment.  However,  in  [11]  the  OVE-ISP  algorithm  did  not  demonstrate  any  volume 
reduction  advantages  over  a  pseudorandom  white  noise  sequence  of  equal  expected 
energy,  or  over  an  alternative  input  synthesis  algorithm  found  in  [40],  when  comparing 
steady-state  performance. 

Therefore,  in  the  time-invariant  case  where  steady-state  performance  is  often  more 
important  than  transient  performance,  other  input  synthesis  strategies  could  be  used. 
However,  in  the  time- varying  case,  where  the  transient  performance  is  very  important, 
OVE-ISP  could  be  used  in  conjunction  with  an  adaptive  robust  controller.  Initially, 
OVE-ISP  could  be  used  to  quickly  reduce  the  ellipsoid  volume  so  that  a  robust  con¬ 
troller  could  be  designed.  It  could  also  be  used  periodically  if  the  ellipsoid  volume 
becomes  to  large  for  effective  robust  control. 

Second,  during  our  course  of  study  of  OVETV  and  OVE-ISP  we  discovered  a 
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significant  property  of  OVE-ISP  which  may  make  it  the  preferred  input  synthesis 
strategy  for  a  certain  class  of  linear  systems  which  may  be  either  time-invariant  or 
time-varying.  Because  of  the  “feedback  nature”  of  OVE-ISP,  the  OVE-ISP  demon¬ 
strates  a  “stabilizing”  property  that  the  pseudorandom  white  noise  or  input  synthesis 
strategy  of  [40]  did  not  possess.  Consequently,  OVE-ISP  is  probably  the  preferred 
strategy  for  systems  which  may  be  unstable  or  marginally  stable. 

Third,  OVE-ISP  offers  design  flexibility  which  may  make  it  a  useful  strategy  not 
only  for  system  identification  but  also  for  EDI  (which  is  discussed  in  the  next  chapter). 
For  an  example  of  FDI-based  input  design,  see  [41]  where  Sadegh,  Madsen,  and  Holst 
discuss  an  optimal  input  design  strategy  for  fault  detection  and  diagnosis  based  on 
stochastic  properties  of  the  system.  In  OVE-ISP,  the  goal  is  to  direct  the  regression 
vector,  such  that  maximum  reduction  in  ellipsoid  volume  occurs.  For  FDI  an 
alternative  goal  might  be  to  direct  (f>k  such  that  fast  tracking  and  isolation  take  place 
after  a  fault  is  detected.  It  should  be  pointed  out,  however,  that  the  OVE-ISP  strategy 
is  highly  dependent  on  the  ellipsoid  center  estimate,  Oko,  which  may  not  adequately 
represent  the  system  dynamics  after  a  fault  occurs. 

Finally,  one  of  the  most  significant  limitations  of  the  OVE-ISP  algorithm  in  [11]  is 
its  lack  of  capability  to  handle  plants  where  the  known  system  delay,  /,  is  greater  than 
1.  This  should  not  be  a  problem  for  physical  systems  where  the  original  continuous 
system  can  adequately  be  modeled  with  rational  functions  in  the  Laplace  domain. 
However,  it  will  not  be  very  satisfactory  for  systems  where  there  is  a  known  trans¬ 
portation  lag  such  that  /  >  1  or  or  /  ^  1.  This  include  systems  such  as  the  crude  oil 
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distillation  column  which  was  discussed  in  Sections  3.4  and  4.2.3.  The  development 
of  OVE-ISP  in  this  section  will  eliminate  this  restriction  and  improve  the  speed  of 
the  OVE-ISP  for  systems  where  I  =  0. 

We  begin  by  developing  the  OVE-ISP  algorithm  for  the  general  case,  /  >  0.  Next, 
we  show  how  the  OVE-ISP  algorithm  can  be  used  in  conjunction  with  OVETV  for 
PSE  of  time- varying  systems.  We  compare  the  performance  of  OVETV  using  OVE- 
ISP  with  the  performance  of  OVETV  using  pseudorandom  white  noise  and  OVETV 
using  the  input  synthesis  strategy  of  [40] .  We  compare  this  performance  for  systems 
which  are  both  stable  and  unstable. 


4.5.1  OVE-ISP  Development  for  General  Delay  Case 

As  discussed  earlier,  the  goal  of  OVE-ISP  is  to  drive  the  system  such  that  <j)k  is 
aligned  with  the  ellipsoid  axis  of  greatest  length.  At  time  this  axis  is  given  by 
the  eigenvector,  ±5^^,  which  corresponds  to  the  largest  magnitude  eigenvalue  of  Pk^ 
and  has  been  normalized  such  that  ||ifc|l2  =  1-  At  some  future  time  ki,  the  desired 
regression  vector,  (j)di  is  selected  to  be  [(pd  —  where  ^  is  a  scalar  which  is 

maximized  subject  to  the  constraints  of  the  physical  system,  and  the  sign  of  Xko  is 
chosen  to  minimize  the  synthesized  input  energy  for  a  given  g. 

Invoking  a  certainty  equivalence  argument,  the  current  ellipsoid  center,  is 
chosen  to  realize  a  state  space  model  of  the  system  which  has  the  form 


<i>k+i  —  4-  Buk-i+i 


(4.64) 
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where 

0,1  •  •  •  Oiji—i  otji  bi  •  •  •  hfYi—i  bfn 

1  0  .  0 

0  •••  i 

;  •••  0  i 

:  1  i 

0  .  0  1  0 

At  this  point  in  our  design,  we  will  assume  that  the  center  ellipsoid  estimate, 
is  equal  to  the  “true”  center  estimate,  that  the  noise  or  uncertainty  entering  our 
system,  Vk,  is  0,  and  that  the  system  is  time-invariant.  Most  likely,  all  three  of  these 
assumptions  will  fail.  However,  the  degree  to  which  they  hold  will  directly  affect  the 
ability  of  OVE-ISP  to  achieve  it  goal. 

Also  at  issue  is  the  observability  and  controllability  of  (4.64).  By  realizing  (4.64) 
as  we  did,  the  states  are  directly  measurable;  therefore  observability  is  not  an  issue. 
If  the  modeling  orders  are  chosen  correctly,  but  (4.64)  is  not  controllable,  it  should 
be  possible  to  choose  another  parameter  estimate  within  Ek^  such  that  (4.64)  is 
controllable.  For  a  more  complete  discussion  of  the  issues  involved  in  choosing  this 
realization,  see  [29]  or  [21]. 

The  system  realization  in  (4.64)  can  now  be  used  to  generate  the  desired  input 
sequence.  Let  ki  =  ko  +  I  —  1.  From  (4.64),  we  know  that 


^ki+i  =  A^ki  +  Buk,-i+i 
=  A(j)ki  +  Buko . 


(4.66) 
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Iterating  (4.66)  and  using  back  substitution,  it  is  easy  to  show  that 
(j^ki+r  =  ^  <t*ki  +  ^  ^BUkQ  +  •  •  •  +  BUkQ.^r-\. 

Equation  (4.67)  can  now  be  rewritten  as 

(4.68) 

where  U  —  Ukg  Uk^+i  •  •  •  Uko+r-i  ]  is  a  vector  which  contains  the  next  r  inputs 
and  (7  =  B  AB  ■  ■  •  A^~^B  is  the  controllability  matrix  which  is  assumed  to 
be  invertible.  Thus,  the  desired  input  sequence  is  given  by 

U  =  C-\<f>d  -  A^(f>k,),  (4.69) 

where  (f>ki+r  has  been  set  to  our  desired  regression  vector,  <j)i. 

Finally,  we  need  to  specify  (f>kt .  It  is  easy  to  show  that 

<l>ko-l  ^  =  0 

<l>ki  -  <  ^ko  ^  =  1  (4.70) 

^  A‘~^(j)ko  +  A^~^Buk()-i+i  +  •  •  •  +  Buko-i  /  >  2, 

where  all  the  terms  on  the  right  hand  side  are  known.  The  desired  input  sequence 

series  takes  r  steps.  The  desired  regressor  is  achieved  at  time  ki  =  ki+r  =  ko+l—i+r. 

This  algorithm  reduces  to  the  algorithm  discussed  in  [11]  for  the  case  when  /  =  1.  It 

is  faster  for  the  case  where  /  =  0.  Finally,  it  handles  the  cases  where  I  >  2  that  are 

not  handled  by  the  development  in  [11]. 

The  steps  of  the  revised  OVE-ISP  algorithm  can  now  be  summarized: 

1.  Use  the  current  OVE  (OVETV)  ellipsoid  center  estimate,  Ok^,  to  generate  the 


state  space  realization  of  (4.64). 
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2.  Find  the  eigenvector  Xk,  which  corresponds  to  the  largest  eigenvalue  of  Fkg.  Set 
(f>j,  =  ±gxka  where  g  is  prespecified. 

3.  Compute  the  input  sequence  specified  by  (4.69)  and  (4.70)  which  is  necessary 
to  transfer  <j)ko  to  <f>d. 

Note  that  even  through  <j)k  will  not  reach  (f>d  until  ki  =  ko  +  I  —  1  +  r,  the  input 
sequence  is  only  specified  until  (fco+»'  —  !)•  Therefore,  if  desired,  a  new  sequence  could 
be  generated  starting  at  {ko  +  r).  If  fact,  if  a  new  sequence  is  to  be  generated  every  r 
steps,  from  (4.69)  we  see  that  the  OVE-ISP  procedure  could  be  viewed  (and  analyzed) 
as  a  multi-rate  time-varying  feedback  scheme.  To  illustrate  this,  see  Figure  26  for 
the  case  where  /  =  1,  and  where  {Ac,Bc,Dc}  represent  the  plant’s  continuous  time 
system  dynamics,  T  is  the  sampling  time,  and  ZOH  is  a  zero-order  hold.  If  an  OVE- 
ISP  sequence  begins  at  time  A;o,  for  k  =  ko  . .  .ko+r-i,  the  time- varying  gain  Kk  is 
given  by  [  •  •  •  ATp+r-i 

Feedback  can  cause  difficulties  for  system  identification  algorithms  [34].  It  is 
shown  in  [42]  how  identifiability  can  be  lost  due  to  static  gain  feedback.  In  this 
case,  the  parameters  cannot  be  uniquely  determined.  However,  identifiability  can  be 
regained  by  using  a  time- varying  feedback  gain  or  feedback  of  sufficiently  high  order. 
The  feedback  gain  in  the  OVE-ISP  procedure,  Kk,  is  time- varying.  Furthermore,  the 
feedback  depends  not  only  on  the  output,  but  on  past  values  of  the  output  and  input. 
While  the  full  implications  of  this  “feedback  nature”  will  require  further  study,  its 
importance  will  be  demonstrated  in  the  next  section. 
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Figure  26:  OVE-ISP  procedure  viewed  as  time-varying  feedback  scheme 


4.5.2  Application  of  OVE-ISP  and  OVETV  to  Time- Varying 
Systems 

As  previously  mentioned,  one  of  the  most  significant  applications  for  OVE-ISP  may 
be  in  the  areas  of  adaptive  PSE  and  adaptive  robust  control.  However,  in  the  devel¬ 
opment  of  the  previous  section,  it  was  explicitly  assumed  that  the  system  we  were 
trying  to  identify  was  time-invariant.  Clearly,  we  would  expect  the  performance  of 
OVE-ISP  to  deteriorate  for  systems  which  vary  “quickly”  with  time.  Consequently, 
for  demonstration  purposes  we  will  work  only  with  “slowly- varying”  systems. 

In  this  section,  we  will  consider  a  variation  of  the  example  examined  in  [40]  and 
[11];  here,  we  will  allow  the  system  parameters  to  vary  with  time.  Consider  the 
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time- varying  system 


where 


and 


t/k  —  ^Ic  4'k  H"  Vh, 
^fc+l  =  ^k  +  '^ki 


0k 

<i>k 


Olfc  (‘■2k  hk  ] 


-yk-1  —yk-2  uk-1 


do  =  [  0.4  0.85  0.75 


Wk 


0.005  -0.002  0.007 


(4.71) 

(4.72) 


(4.73) 

(4.74) 


(4.75) 

(4.76) 


and  Vk  is  a  uniformly  distributed  random  variable  between  [-0.05  0.05].  The  system 
is  initially  at  rest. 

To  apply  OVETV,  the  bounds  in  (2.3)  and  (2.4)  must  be  specified.  We  set  Tk  = 
— =  0.05  and  Rk  =  where  Cfc  =  ||infc||2  and  I  is  the  identity  matrix.  The 
orientation  matrix  Pq  is  set  to  10^7,  and  the  center  9q  is  set  to  the  origin. 

For  this  example,  we  will  use  three  different  input  sequences.  OVE-ISP  with  q  =  \ 
is  repeatedly  applied  to  the  system  to  generate  the  first  input  sequence.  The  second 
input  sequence  is  a  uniformly  distributed  random  white  noise  sequence  whose  bounds 
are  set  such  that  it  has  the  same  input  energy  as  the  OVE-ISP  sequence.  The  final 
input  sequence  uses  an  approach  which  was  developed  by  Pronzato  and  Walter  and 
can  be  found  in  [40].  This  input  sequence  is  “bang-bang”  in  nature  and  is  designed 
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Figure  27:  First  time- varying  example  comparing  system  inputs 


using  an  “optimal”  one-step-ahead  cost  function  with  a  constraint  on  the  magnitude 
of  the  input.  We  will  choose  the  bound  on  the  input  level  so  that  the  energy  of  this 
input  sequence  also  matches  the  energy  of  the  OVE-ISP  sequence. 

Before  applying  any  of  these  input  sequences,  a  uniformly  distributed  random 
white  noise  sequence  of  length  3  was  applied  to  the  system  to  “initialize”  the  OVETV 
algorithm.  The  ellipsoid  volumes  for  the  three  input  sequences  are  shown  in  Figure  27. 
While  OVE-ISP  may  have  a  slight  advantage  during  the  transient  phase,  overall,  there 
appears  to  be  no  clear  advantage  for  any  of  the  input  sequences.  The  parameters, 
center  estimates,  and  upper  and  lower  bounds  are  shown  in  Figure  28  for  the  OVE-ISP 
input  sequence.  The  plots  for  the  other  input  sequences  were  similar. 
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Figure  28;  First  time-varying  example  with  OVE-ISP  input  sequence 

This  simulation  is  repeated  as  above  except  with  Wk  now  given  as 
Wk  =  [  0.005  0.0035  —0.007  .  The  ellipsoid  volumes  for  this  second  example 
are  shown  in  Figure  29.  Upon  first  inspection,  it  would  appear  that  the  performance 
of  OVE-ISP  is  inferior  to  the  performance  of  the  random  input  and  the  Pronzato- 
Walter  algorithm.  However,  if  we  examine  the  output  sequences  in  Figure  30  for  each 
of  the  three  input  input  sequences,  we  clearly  see  the  source  of  the  discrepancy.  The 
output  sequences  for  the  random  input  and  Pronzato- Walter  algorithms  are  growing 
to  levels  which  would  probably  be  unacceptable  for  most  physical  systems,  while  the 
output  sequence  for  the  OVE-ISP  algorithm  is  remaining  “well-behaved.”  Clearly,  the 
energy  of  the  output  sequence  has  as  much  effect  on  the  volume  of  the  ellipsoids  as 


4.5.  Input  Synthesis 


103 


Figure  29:  Second  time- varying  example  comparing  system  inputs 


the  energy  of  the  input  sequence. 

The  reason  for  the  unacceptable  increase  in  magnitude  of  the  output  sequences 
for  the  random  input  and  Pronzato- Walter  algorithm  can  be  seen  in  Figure  31.  In 
Figure  31,  we  see  that  while  the  system  poles  at  each  time-step  during  the  first 
example  remain  within  the  unit  circle,  the  system  poles  during  the  second  example 
clearly  move  outside  the  unit  circle.  While  the  locations  of  the  “frozen-time”  poles  by 
themselves  do  not  guarantee  stability  or  instability  for  linear  time- varying  systems, 
they  do  provide  information,  particularly  for  slow-varying  systems,  about  the  local 
(in-time)  behavior  of  the  system. 

In  the  second  example,  once  the  poles  of  the  system  move  outside  the  unit  circle. 
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Figure  30:  System  outputs  for  second  example 


Figure  31:  Frozen-time  poles  of  first  and  second  examples 
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Figure  32:  Second  time-varying  example  with  OVE-ISP  input  sequence 


the  magnitude  of  the  output  of  the  open-loop  system  grows  rapidly.  However,  the 
output  sequence  using  the  OVE-ISP  input  sequence  remains  well-behaved  because  of 
the  “feedback  nature”  of  its  synthesis.  Furthermore,  as  seen  in  Figure  32,  OVETV  is 
able  to  track  the  parameters  effectively  using  the  OVE-ISP  input  sequence.  This  is 
in  contrast  to  the  Pronzato- Walter  case  where  the  overall  volume  is  smaller,  but  the 
uncertainty  around  the  parameter  bi  becomes  very  large,  as  can  be  seen  in  Figure  33. 
OVE-ISP  has  achieved  similar  performance  for  unstable  linear  time-invariant  systems. 
In  conclusion,  while  this  “stabilizing”  property  is  appealing,  clearly  more  analysis  is 
needed  to  determine  when  and  if  this  property  can  be  guaranteed. 
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Figure  33:  Second  time-varying  example  with  Pronzato- Walter  input  sequence 

4.6  Summary 

In  this  chapter,  we  discussed  several  algorithm  modifications  and  extensions  which 
serve  to  improve  one  or  more  properties  of  the  OVE  or  OVETV  algorithms.  In  Sec¬ 
tion  4.2,  we  showed  how  MISO  systems  could  be  placed  in  the  OVE  or  OVETV  frame¬ 
work.  We  began  by  showing  how  linear  time-invariant  MISO  state  space  equations 
can  be  placed  in  the  OVE  framework.  Next,  we  extended  Theorem  5  to  linear  time- 
varying  MISO  systems  which  satisfy  a  certain  “observability”  property.  Theorem  6 
showed  how  systems  which  satisfy  this  property  can  be  transformed  into  the  OVETV 
framework.  Finally,  results  were  presented  showing  how  the  OVETV  algorithm  could 
be  applied  to  the  crude  oil  distillation  column  using  the  MISO  framework. 
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In  section  4.3,  we  showed  how  dependencies  in  parameter  variations  could  be  ex¬ 
ploited  to  reduce  the  number  of  calculations  necessary  to  implement  the  OVETV 
algorithm.  When  these  dependencies  are  present,  support  functions  are  utilized  to 
handle  the  degeneracies  which  arise  in  the  ellipsoid  bound  on  the  parameter  distur¬ 
bance  vector,  w/..  Properties  of  determinants  are  used  to  decrease  the  size  of  the 
generalized  eigenvalue  problem  in  step  1  of  the  optimal  time  update  equations,  as 
well  as  the  size  of  the  polynomial  which  needs  to  be  solved  in  step  2. 

Section  4.4  introduced  a  square-root  implementation  of  the  OVE  algorithm.  This 
implementation  serves  to  prevent  the  ellipsoid  orientation  matrix,  Pk,  from  becoming 
non-positive  definite  due  to  finite  precision  calculations.  The  implementation  fol¬ 
lows  an  outline  similar  to  Potter’s  algorithm  [36]  which  is  often  used  in  stochastic 
estimation  routines.  A  comparison  with  the  direct  implementation  of  the  OVE  algo¬ 
rithm  showed  that  while  the  square-root  implementation  requires  more  floating  point 
operations,  the  number  of  operations  is  of  the  same  order  for  each  implementation. 
Furthermore,  in  most  cases  this  increase  in  number  of  computations  will  be  less  im¬ 
portant  than  the  increase  in  numerical  stability  which  is  gained  by  the  square-root 
implementation.  For  the  time- varying  case  it  was  shown  how  the  scalar  addition  and 
scalar  multiplication  algorithms  of  Section  2.6  can  also  be  implemented  in  a  square- 
root  framework. 

Finally,  in  section  4.5  the  OVE-ISP  input  synthesis  procedure  of  [21]  and  [11] 
was  investigated.  First,  we  extended  OVE-ISP  to  handle  the  class  of  systems  which 
contains  a  known  transportation  lag.  The  modified  algorithm  also  improves  the  speed 
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of  response  for  the  case  where  there  is  no  delay.  Next,  we  demonstrated  how  OVE-ISP 
and  OVETV  could  be  applied  to  linear  time- varying  systems.  The  simulation  results 
presented  demonstrate  a  “stabilizing”  property  of  OVE-ISP  which  is  not  possessed  by 
a  random  input  sequence  and  an  alternative  synthesis  procedure  found  in  [40].  This 
“stabilizing”  property  makes  OVE-ISP  the  preferred  input  scheme  for  unstable  linear 
time-invariant  and  linear  time- varying  systems. 


CHAPTER  V 


Fault  Detection  and  Isolation  using  PSE 

5.1  Overview 

To  attain  autonomy  in  control  systems  it  is  necessary  that  a  controller  be  able  to 
detect  and  identify  faults  in  a  complex  dynamical  system.  These  faults  (or  failures) 
may  occur  in  the  sensors,  actuators,  and  components  of  the  system  we  are  trying  to 
control.  It  is  very  important  that  a  fault  detection  and  isolation  (FDI)  scheme  be  able 
to  accurately  detect  failures  in  the  presence  of  modeling  errors,  system  and  measure¬ 
ment  noise,  and  parameter  variations.  Because  of  these  stringent  requirements,  this  is 
a  natural  setting  for  a  robust  system  identification  technique  such  as  set-membership 
identification.  In  this  chapter  we  show  how  the  OVE  the  OVETV  algorithms  can  be 
used  for  detection  and  isolation  of  faults  in  dynamical  systems. 

The  chapter  begins  with  a  brief  overview  of  FDI  and  some  of  the  issues  involved 
in  its  implementation.  Next,  we  will  discuss  two  methods  for  detecting  faults  in 
dynamical  systems.  The  first  method  relies  on  the  consistency  check  which  is  integral 
to  the  OVE  and  OVETV  algorithms,  while  the  second  method  utilizes  an  ellipsoid 
intersection  test  to  detect  a  fault. 
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If  a  fault  is  detected  by  an  inconsistency  check  in  the  OVE  or  OVETV  algorithms, 
a  fault  is  signalled  and  the  algorithms  halt.  However,  in  many  situations  it  is  desir¬ 
able  to  track  the  parameters  in  the  parameter  space  after  a  fault  is  indicated.  Two 
“recovery”  strategies  are  presented.  The  first  strategy  simply  resets  the  current  ellip¬ 
soid  to  a  “large  enough”  ellipsoid  guaranteed  to  contain  the  “true”  parameter  after  a 
fault  is  detected.  This  is  the  only  method  that  we  know  of  that  can  guarantee  that 
the  “true”  parameter  will  be  contained  in  the  parameter  set  immediately  after  a  fault 
it  indicated.  However,  after  the  ellipsoid  is  initially  reset,  there  may  be  a  transition 
time  for  which  the  ellipsoid  is  too  large  to  be  effective  for  applications  such  as  robust 
adaptive  control. 

Consequently,  an  alternative  strategy  is  also  introduced.  This  strategy  uses  a 
projection  strategy  to  combine  the  information  contained  in  the  new  measurements 
with  the  “best”  of  the  information  retained  by  the  previous  ellipsoid.  While  this 
strategy  is  not  guaranteed  to  capture  the  “true”  parameter,  it  often  succeeds  in 
doing  so  very  quickly.  Rules  are  also  proposed  for  combining  these  strategies  to  take 
advantage  of  the  guaranteed  properties  of  the  resetting  algorithm  and  the  transition 
properties  of  the  projection  algorithm. 

Finally,  the  algorithms  of  this  chapter  are  illustrated  using  the  linear  time- varying 
circuit  example  of  Section  3.3.  Some  techniques  for  fault  isolation  are  also  discussed 
in  our  concluding  remarks. 
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5.2  Issues  in  FDI 

The  use  of  FDI  is  becoming  more  and  more  critical  as  the  complexity  of  systems  to 
be  controlled  increases.  Applications  where  FDI  is  necessary  include  complex  process 
plants  where  early  detection  of  slowly- varying  faults  can  avoid  major  plant  break¬ 
downs  and  disasters.  FDI  is  also  applicable  to  high  performance  vehicles,  such  as 
ships,  submarines,  airplanes,  and  spacecraft,  where  safety  and  significant  financial  in¬ 
vestment  are  at  stake.  The  application  of  FDI  techniques  has  also  been  made  possible 
by  the  increased  computational  capability  of  digital  computers  and  the  development 
of  advanced  processing  techniques. 

There  are  four  basic  approaches  to  FDI:  (1)  limit  and  trend  checking,  (2)  physical 
redundancy,  (3)  analytical  redundancy,  and  (4)  knowledge-based  redundancy.  Limit 
and  trend  checking  is  the  simplest  and  the  most  widely  used.  In  this  strategy,  limits 
are  placed  on  measured  variables  or  on  there  rate  of  change.  If  these  limits  are 
surpassed,  a  fault  is  detected.  Examples  include  smoke  alarms  and  “dummy”  lights 
on  an  automobile. 

With  physical  redundancy,  multiple  sensors  are  used  for  each  measurement.  A 
voting  strategy  is  used  to  decide  which  measurements  to  use.  Multiple  actuators  may 
also  be  used.  While  this  method  has  the  advantage  of  being  simple,  the  additional 
hardware  requirements  may  be  costly  in  terms  of  dollars  as  well  as  mass  and  volume 
requirements. 

Analytical  redundancy  is  an  approach,  often  model-based,  which  generates  residu¬ 
als  utilized  to  detect  and  isolate  failures.  This  is  the  approach  used  in  this  dissertation 
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and  will  be  discussed  in  detail.  Knowledge-based  redundancy  uses  heuristic  model 
information  and  often  functions  as  a  supervisor  utilizing  the  other  approaches. 

One  model  considered  by  analytical  redundancy  techniques  is 

x{{k  +  l)T)  =  Ax{kT)  +  Bu{kT)  +  Ed{kT)  +  Kf{kT)  (5.1) 

y{kT)  =  Cx{kT)  +  Fd{kT)  +  Gf{kT),  (5.2) 

where  a;  G  is  the  state,  «  €  is  the  control  input,  j/  €  in  the  output,  d  G  3?^ 
is  the  unknown  input  vector,  /  G  St*”  is  the  fault  vector,  and  T  is  the  sampling  time. 
The  unknown  input  vector,  d,  represents  noise  and  modeling  errors  which  may  cause 
false  alarms  in  a  FDI  strategy.  The  fault  vector,  /,  is  nonzero  only  when  a  fault 
occurs. 

Another  type  of  model  considered  by  analytical  redundancy  techniques  is 

yk  =  0l<f>k  +  Vk,  (5.3) 

where  yk  is  the  system  output,  9k  is  the  parameter  vector,  (f)k  is  the  regression  vector, 
and  Vk  is  the  system  disturbance.  Faults  are  modeled  by  changes  in  the  parameter 
vector,  9k,  which  are  not  anticipated  by  the  no-fault  behavior  of  the  system. 

There  are  two  types  of  fault  modes  with  which  we  are  concerned;  abrupt  and 
incipient.  Abrupt  faults  are  step-like  changes  which  require  quick  detection  for  safety 
reasons.  Incipient  faults  are  slowly  varying  failures,  e.g.,  bias  or  drift,  which  take 
longer  to  detect,  but  are  usually  more  maintenance  related. 

As  discussed  previously,  analytical  redundancy  relies  on  the  generation  of  residuals 
to  detect  and  isolate  faults.  Residuals  are  simply  functions  which  are  accentuated  by 
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the  fault  vector  /  and  are  ideally  zero  when  f  —  0.  Most  analytical  redundant 
schemes  require  two  steps:  (1)  generation  of  residuals  and  (2)  decision  and  isolation 
of  the  faults.  Residuals  are  most  often  generated  using  either  state  estimation  or 
parameter  estimation.  Examples  of  state  estimation  schemes  include  parity  checks, 
observer  schemes,  and  detection  filters.  Once  residuals  are  generated,  they  may  be 
used  to  form  decision  functions  where  the  relevant  information  is  more  apparent. 
Isolation  of  the  failures  is  based  on  a  fault  signature,  which  is  a  signal  defining  the 
effects  associated  with  a  particular  fault  (often  derived  from  a  model  of  the  faulty 
system). 

To  achieve  fault  detection  and  isolation,  three  types  of  models  may  be  used:  nom¬ 
inal  representing  no  fault  behavior,  actual  representing  observed  behavior,  and  faulty 
representing  the  system  after  a  fault  has  occurred.  Ultimately,  we  may  desire  to  know 
the  type,  size  and  source  as  well  as  the  time  and  location  of  the  fault.  To  achieve  this, 
often  additional  heuristic  information  is  required,  i.e.,  knowledge  based  redundancy. 

While  FDI  has  received  much  attention  in  the  research  community  (see  e.g.,  [43], 
[44],  and  [45]),  there  has  been  little  investigation  into  the  use  of  set-membership  iden¬ 
tification  approaches  for  FDI.  This  is  surprising  given  the  stringent  requirements  for 
robustness.  Set-membership  identification  algorithms  have  the  very  desirable  prop¬ 
erty  that  the  estimated  parameter  sets  are  guaranteed  to  contain  the  “true”  param¬ 
eter.  It  should  be  mentioned,  though,  that  this  robustness  property  is  only  as  good 
as  the  assumptions  upon  which  the  parameter  estimation  algorithms  are  based. 

Work  which  is  related  to  the  use  of  PSE  for  FDI  can  be  found  in  [46],  where  a 
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two-ellipsoid  overlap  test  for  on-line  failure  detection  is  developed.  The  authors  did 
not  use  set-membership  identification  but,  instead,  used  confidence  regions  based  on 
covariance  matrices.  They  compared  the  confidence  region  of  an  on-line  identified 
model  of  the  system  with  the  confidence  region  of  a  nominal  model.  A  fault  was 
detected  when  the  confidence  regions  did  not  overlap. 

In  [14],  the  authors  show  how  a  version  of  the  OBE  algorithm  developed  in  [4]  can 
be  modified  by  artificially  increasing  the  noise  bound  to  guarantee  consistency  if  an 
update  takes  place.  This  algorithm  suffers  from  the  fact  that  this  consistency  can  still 
be  lost  if  no  update  takes  place.  However,  they  have  developed  a  “rescue  procedure” 
for  when  tracking  is  lost,  for  example  due  to  an  unexpected  large  parameter  jump 
which  could  be  caused  by  a  fault.  Other  rescue  procedures  are  discussed  in  [47]. 

5.3  Detection  of  Failures 

A  critical  part  of  an  EDI  scheme  is  the  ability  to  detect  failures.  In  this  section, 
we  will  discuss  two  approaches  for  detection  of  failures.  The  first  approach  uses 
the  consistency  check  which  is  integral  to  the  OVE  and  OVETV  algorithms.  This 
check  determines  if  the  intersection  between  the  ellipsoid  resulting  from  the  time 
update,  Ek,  and  the  region  based  on  the  new  measurements,  Fk,  is  empty.  The  second 
approach  combines  the  ellipsoid  intersection  test  of  [46]  with  the  OVETV  algorithm. 
An  alternative,  more  straight-forward,  development  of  the  ellipsoid  intersection  test 
is  given. 
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5.3.1  Consistency  Check 

Consider  the  system  given  by  (2.1)  and  (2.2)  with  known  bounds  given  in  (2.3)  and 
(2.4).  To  apply  the  OVETV  algorithm  to  this  system,  we  assume  that  the  following 
are  known:  the  system  structure,  i.e.,  n,  m,  and  /,  the  bounds  on  the  system  distur¬ 
bance,  i.e.,  7^  and  7^,  the  bound  on  the  parameter  disturbance  vector,  i.e.,  Rk,  and 
the  initial  ellipsoid  which  is  guaranteed  to  contain  the  “true”  parameter  vector,  i.e. 
Eq.  Under  the  assumption  that  these  parameters  are  known,  the  algorithm  produces 
an  ellipsoid  at  each  time  k,  Ek,  which  is  guaranteed  to  contain  the  “true”  parameter 
vector. 

However,  if  any  of  these  are  unknown,  we  may  have  the  case  where  EkCl  Fk  =  0, 
where  Fk  is  the  region  between  the  two  parallel  hyperplanes,  and  Ek  is  the  ellipsoid 
resulting  after  the  time  update  equations.  For  an  example,  consider  the  first-order 
system  in  Section  3.2.  With  the  algorithm  initialized  as  before,  at  A;  =  16  a  fault  is 
induced  causing  the  parameter  Uk  to  jump  by  0.25,  i.e.,  aie  =  cti5 -1-0.25.  This  violates 
the  bounds  on  the  parameter  disturbance  vector,  Wk,  and  produces  the  simulation 
results  shown  in  Figure  34.  Clearly,  Eie  Cl  Fie  =  0,  indicating  that  a  fault  has 
occurred. 

Let  us  formalize  an  approach  implied  by  this  simple  example.  Let  the  nominal 
system  model  be  given  by 


Vk  —  &k^k  +  Vkt 

(5.4) 

Ok+l  ^  Ok  +  Wk, 

(5.5) 
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Figure  34:  Fault  detected  at  A:  =  16 


where  6k 


<j>k  = 


-Vk-i 


aik  • •  • 

Vk-n  Uk-l 


^nk  6lk  ■  ■  ‘  bmk 

lT 


,  6q  €  Eq^  and 
.  The  system  disturbance,  Ufc,  and 


the  parameter  disturbance  vector,  Wk,  are  assumed  to  satisfy  the  following  bounds 


Xk<Vk<  lki  (5-6) 

{wk  €  r :  e'^wk  <  \JeTRk6-,  6  €  r },  (5.7) 

where  7^  <  7^^,,  and  ^  €  9?,  7*,  €  and  Rk  €  are  known  at  each  time  k.  The 
matrix  Rk  is  symmetric,  semi-positive  definite. 


From  these  quantities,  define  the  following  sequences: 

Un  =  {uk}k=0,...,N  Yn  =  {yk}k=0,...,N 

Vn  =  {vk}k=0,...,N  Wn  =  {wk}k=0,...,N- 


(5.8) 
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For  this  work,  we  will  define  a  fault  as  follows. 

Definition  1  For  a  given  input  sequence,  Un,  o-nd  output  sequence,  Y^,  of  the  actual 
system,  a  fault  is  said  to  have  occurred  if  there  does  not  exist  a  pair  of  sequences,  Vn 
and  Wn,  which  satisfy  the  nominal  system  equations,  (5.4)  and  (5.5),  and  nominal 
system  bounds,  (5.6)  and  (5.7). 

We  can  now  state  the  following  theorem. 

Theorem  7  Given  that 

•  the  actual  system  satisfies  (5.4)-(5.7)  under  nominal  operating  conditions, 

•  the  OVETV  algorithm  is  initialized  such  that  7jt  >  7^,  7^  <  Rk  >  Rk,  and 
Eo  D  Eo,  and 

•  the  OVETV  algorithm  is  applied  to  the  input  sequence,  Un,  and  output  sequence, 
Yn,  of  the  actual  system. 

If 

Ekr\Fk  =  0,  (5.9) 

for  any  k  E  {0, . . . ,  N},  then  a  fault  has  occurred. 

Proof:  Given  the  input  sequence,  and  output  sequence,  YJv,  of  the  actual  system, 
assume  that  Ek  Fk  =  0,  but  that  no  fault  has  occurred.  Therefore,  there  exist 
sequences,  Vn  and  Wn,  such  that  the  nominal  system  equations  and  nominal  systems 
bounds  are  satisfied.  However,  if  Vn  and  Wn  exist  which  satisfy  (5.4)-(5.7),  by 
construction  of  the  OVETV  algorithm. 
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A 

This  is  a  contradiction.  Therefore,  if  EkCFk  =  0,  then  there  cannot  exist  sequences, 
Vn  and  Wn,  such  that  the  nominal  system  and  nominal  system  bounds  are  satisfied. 
Thus,  a  fault  must  have  occurred.  □ 

In  this  section  we  applied  Theorem  7  to  the  first-order  example  of  Section  3.2. 
In  Section  5.5,  we  will  demonstrate  this  detection  strategy  on  a  more  complicated 
example,  the  linear  time-varying  circuit  of  Section  3.3. 

When  a  fault  is  detected  by  this  algorithm,  i.e,  when  (1  =  0,  the  OVETV 

algorithm  halts.  However,  there  are  many  situations  where  it  would  be  desirable  to 
track  the  parameters  after  a  fault  is  detected.  In  Section  5.4,  algorithms  are  discussed 
which  allow  the  OVETV  algorithm  to  “recover”  after  a  fault  is  detected.  In  the  next 
section,  we  will  discuss  another  algorithm  for  fault  detection.  This  algorithm  has  the 
advantage  that  for  linear  time-invariant  systems  with  incipient  faults,  no  “recovery” 
strategy  is  needed  to  continue  tracking  the  parameters  after  a  fault  is  detected. 

5.3.2  Ellipsoid  Intersection  Test 

Consider  the  system  given  by  (2.1)  and  (2.2)  with  known  bounds  given  in  (2.3)  and 
(2.4).  Assume  that  under  nominal  conditions,  the  system  can  be  described  by  (5.4)- 
(5.7)  where  Rk  is  known  to  be  equal  to  0  for  all  k,  i.e.,  the  nominal  system  is  time- 
invariant.  If  we  apply  the  OVE  algorithm  to  the  system  when  it  is  operating  under 
nominal  conditions,  we  will  get  a  resulting  ellipsoid,  Enom- 

To  monitor  the  system  for  incipient  (slowly-varying)  faults,  we  can  apply  the 
OVETV  algorithm  with  Rk  set  greater  than  0  such  that  parameter  tracking  is  never 
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lost.  In  this  situation,  we  can  compare  the  ellipsoid,  Ek,  produced  by  the  OVETV 
algorithm  during  actual  operation  with  the  nominal  ellipsoid,  Enom,  which  was  pro¬ 
duced  during  nominal  operating  conditions.  If  these  ellipsoids  intersect,  it  is  possible 
that  the  actual  plant  matches  the  nominal  plant  model.  However,  if  the  ellipsoids  do 
not  intersect,  clearly  a  fault  has  occurred. 

For  an  example,  consider  the  first-order  system  in  section  3.2.  Initially,  the  OVE 
algorithm  was  applied  to  the  system  under  nominal  conditions,  i.e.,  Uk  =  0.5  and 
bk  =  1.0.  This  resulted  in  the  ellipsoid 


Enom  =  {e:{e-  0nomfp-,lie  -  Onom)  <  F,  B  E  ^ (5.11) 


where 


0.0007  0.0006 
0.0006  0.0011 


0.4970 

1.0045 


(5.12) 


During  the  actual  operation  of  the  plant,  the  parameters  begin  to  drift,  just  as  in 
Section  3.2.  In  order  to  track  the  parameters  as  they  drift,  the  OVETV  algorithm 
was  initialized  the  same  as  in  Section  3.2.  In  Figure  35,  we  see  at  that  time  k  =  99, 
the  nominal  ellipsoid,  Enom,  and  the  actual  ellipsoid,  Ek,  fail  to  intersect.  Clearly,  a 
fault  has  occurred. 

In  this  section,  we  will  develop  an  algorithm  to  automatically  detect  if  two  ellip¬ 
soids  intersect.  While  such  an  algorithm  was  previously  developed  in  [46],  the  devel¬ 
opment  in  this  section  serves  to  illustrate  the  process  and  is,  we  feel,  more  straight 
forward  than  that  of  [46].  Furthermore,  we  will  recommend  an  additional  step  to  the 
algorithm  presented  in  [46]  which  often  results  in  a  reduced  number  of  computations. 
Finally,  the  detection  strategy  discussed  above  will  be  stated  in  terms  of  a  theorem. 
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a 


Figure  35:  Fault  detected  at  A:  =  99 


Ellipsoid  Intersection  Equations 

We  wish  to  develop  a  test  to  determine  if  two  ellipsoids, 

El  =  {e:{e-  difPT^e  -  Oi)  <  i;  e  (5.13) 

E2  =  {e:{e-e2fP2\o-e2)<i]02eW}  (5.14) 

intersect.  This  is  equivalent  to  testing  whether  the  following  ellipsoids  intersect 

El  =  {x:x'^Pi^x<l-,xe^"}  (5.15) 

E2  =  {x  :  {x  —  X2)^P2^{x  —  X2)  <1]X2  (5.16) 
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Figure  36:  Intersecting  and  nonintersecting  ellipsoids 


where  x  =  O  —  Oi  and  X2  —  O2  —  O1.  A  pair  of  intersecting  and  a  pair  of  nonintersecting 
ellipsoids  are  shown  in  Figure  36. 

If  the  center  of  Ei  is  contained  within  £2, 


xlP2^X2  <  1,  (5.17) 

then  the  ellipsoids  intersect,  i.e.,  E\C\E2  ^  0.  If  E\r\E2  ^  0  and  (5.17)  is  not  satisfied, 
then  clearly  there  must  be  at  least  one  point,  a;*,  which  lies  on  the  hypersurface  of 
E2,  i.e.,  {x*  —  X2)^P2^{x*  —  X2)  =  1,  and  lies  within  Ei,  i.e.,  x*^P{^x*  <  1. 

This  suggests  the  following  development,  as  a  straight-forward  alternative  to  that 
of  [46].  If  (5.17)  is  satisfied,  then  the  ellipsoids  intersect  and  the  algorithm  stops.  If 


122 


CHAPTER  V.  Fault  Detection  and  Isolation  using  PSE 


(5.17)  is  not  satisfied,  minimize 

J{x)  =  (5.18) 

subject  to  the  constraint  that  x  satisfies 

{x  —  X2)^P2^{x  —  X2)  =  1.  (5.19) 

If  J{x*)  <  1,  where  x*  is  the  optimum  value,  then  the  ellipsoids  intersect,  otherwise, 
they  do  not.  The  value  of  x*  can  be  seen  in  Figure  36  for  a  pair  of  intersecting 
ellipsoids  and  a  pair  of  nonintersecting  ellipsoids. 

To  solve  (5.18)  subject  to  (5.19),  we  form  the  augmented  cost  function 

Ja{x,  A)  =  x^Pi^x  +  A((a;  —  X2)^P2^{x  —  X2)  —  1),  (5.20) 

where  A  is  the  Lagrange  multiplier.  Taking  the  gradient  of  (5.20)  with  respect  to  x 
and  setting  it  equal  to  0,  we  get 

^{x*,X*)  =  2PE^x*  +  2X*Pp\x*  -  X2)  =  0.  (5.21) 

ox 

Rearranging  (5.21)  results  in 

{P{^  +  A*P2'^)x*  =  X*P^^X2.  (5.22) 

Assuming  that  A*  is  positive,  then  the  term  multiplying  x*  is  positive  definite  and 
invertible,  therefore, 

X*  =  X*{Pp^  +  X*P2^)-^P2^X2.  (5.23) 

Taking  the  gradient  of  (5.20)  with  respect  to  A  and  setting  it  equal  to  0  gives 


^(X*,  A*)  =  (X*  -  X2)^P2  -  X2)  -  1  =  0. 


(5.24) 
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By  substituting  (5.23)  into  (5.24)  and  performing  some  algebra,  we  get 


+  A*P2"')"'A"'^2  =  1.  (5.25) 


To  show  that  a  unique  positive  solution  exists  for  A*,  we  make  the  following 
substitution.  Because  Pf^  and  P^^  are  symmetric,  positive  definite,  there  exists  a 
nonsingular  matrix  S  such  that 


Pf^  =  s^s 

(5.26) 

PF^  =  S'^DS, 

(5.27) 

where  D  =  diag(di,  ^2, . . . ,  d^)  [32].  Substituting  (5.26)  and  (5.27)  into  (5.25)  results 


in 


x^S^il  +  X*D)-^D{I  +  X*D)-^Sx2  =  1 


(5.28) 


The  left  hand  side  of  (5.28)  can  be  written  as  rj{X*)  =  ^  tj — ^  where  v  = 

i=i  (1  +  A*j 

Ui  V2  •••  Ur  ]•  Clearly,  r){oo)  =  0.  At  A*  =  0,  7/(0)  =  X2PF^X2  must  be 


SX2  = 


greater  than  1,  otherwise,  we  would  have  stopped  the  algorithm  earlier  because  we 
would  have  known  that  the  ellipsoids  intersect.  For  A*  >  0,  ^(A*)  <  0  and,  therefore, 
?/(A*)  is  strictly  decreasing.  Consequently,  a  unique  positive  solution  must  exist  for 
A*. 

Therefore,  the  algorithm  for  detecting  when  two  ellipsoids,  Pi  and  P2>  intersect 
is  given  as  follows: 


1.  If  X2P2^X2  <  1  where  X2  =  O2  -  6i,  then  Pi  fl  P2  ^  0  and  the  algorithm  stops; 


otherwise  continue. 
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2.  Solve  /(A)  =  xJP2^X3  —  1  for  the  unique  positive  solution,  A*,  where  Xs  = 

{PT^  +  Pp^ {O2  —  Oi).  The  authors  in  [46]  recommend  using  the  bisection 

method  to  solve  for  A*.  This  method  is  guaranteed  to  converge  if  the  initial 
values  for  A  are  set  at  0  and  some  positive  value  such  that  /(A)  <  0. 

3.  If  x*^P2^x*  <  1  where  x*  =  A*(Pf^  +  X*  P2^)~^  P^^  {62  —  Oi)  and  A*  is  found 
as  above,  then  Ei  C]  E2  ^  0,  otherwise,  Ei  0  E2  =  0,  i.e.,  Ei  and  E2  do  not 
intersect. 

We  also  recommend  an  additional  step,  which  is  not  included  in  [46],  to  precede 
step  1. 

0.  If  x^PT^X2  <  1  where  X2  =  O2  —  then  EiC\  E2^%  and  the  algorithm  stops. 

Step  0  determines  if  the  center  of  E2  is  contained  in  E\.  While  this  step  is  not 
necessary  for  the  algorithm  to  achieve  its  goal,  it  will  reduce  the  number  of  times  the 
search  in  step  2  has  to  be  conducted. 


Ellipsoid  Intersection  Detection  Strategy 

Now,  let  us  formalize  the  detection  strategy  which  utilizes  the  ellipsoid  intersection 
test  developed  above.  Let  the  nominal  system  model  be  given  by  (5.4)-(5.7)  with 
Rk  =  0.  Let  the  actual  system  under  all  conditions  be  given  by  (5.4)-(5.6)  and 

{wk  e  r :  e'^wk  <  ^JeTRi^e-  0  e  r },  (5.29) 
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where  Rk  €  is  known  at  each  time  k.  The  matrix,  is  symmetric,  semi-positive 
definite.  We  can  now  state  the  following  theorem. 


Theorem  8  Given  that 

•  the  actual  system  satisfies  (5.4)-(5.7)  under  nominal  operating  conditions  with 

Rk  —  0; 

•  an  ellipsoid,  Enom,  was  produced  by  the  OVE  algorithm  which  was  initialized 
such  that  7^  >  7*.,  7^  <  7^,  and  Eq  3  Eq,  and  which  was  applied  to  the  system 
under  nominal  operating  conditions; 

•  the  OVETV  algorithm  is  initialized  such  that  7^.  >  7^.,  7^  <  7^,  Rk  >  Rk)  o,f^d 
Eq  D  Eq;  and, 

•  the  OVETV  algorithm  is  applied  to  the  input  sequence,  Un,  (‘■nd  output  sequence, 
Yn,  of  the  actual  system. 

If 

EknEnom  =  ^,  (5.30) 

for  any  k  G  {0, . . . ,  N},  then  a  fault  has  occurred. 

Proof:  Given  the  input  sequence,  ^'Hd  output  sequence,  Yn,  of  the  actual  system, 
assume  that  Ek  H  Enom  =  0,  hut  that  no  fault  has  occurred.  Therefore,  there  does 
exist  sequences,  Vjv  and  Wn^  such  that  the  nominal  system  equations  and  nominal 
systems  bounds  are  satisfied.  Therefore,  the  true  parameter,  6*,  must  be  in  Ek  for 
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all  k  —  0, . . .  ,N.  By  construction  of  the  OVE  algorithm,  6*  must  also  be  contained 
within  Enom-  Therefore, 

Ek  n  Enom  ^  0.  (5.31) 

This  is  a  contradiction.  Therefore,  if  EkCEnom  =  0?  then  there  cannot  exist  sequences, 
Vn  and  Wn,  such  that  the  nominal  system  and  nominal  system  bounds  are  satisfied. 
Thus,  a  fault  must  have  occurred.  □ 

Clearly,  this  algorithm  has  an  advantage  over  the  consistency  check  detection 
strategy  in  the  fact  that  for  incipient  faults,  no  recovery  strategy  is  needed.  However, 
it  does  require  calculations  in  addition  to  those  of  the  OVETV  algorithm.  Further¬ 
more,  it  would  be  difficult  to  apply  this  algorithm  to  systems  which  are  nominally 
time- varying,  because  it  would  be  required  to  know  Enom  f’S  a  function  of  k  for  all  time. 
Finally,  abrupt  failures  will  still,  most  likely,  cause  an  inconsistency  in  the  OVETV 
algorithm.  Consequently,  for  abrupt  failures,  no  matter  which  detection  strategy  we 
use,  the  algorithm  recovery  strategies  of  the  next  section  will  be  important. 

5.4  Algorithm  Recovery  Strategies 

When  EkC\  Fk  =  0,  the  OVE  and  OVETV  algorithms  indicate  an  inconsistency  and 
stop.  This  is  satisfactory  if  the  system  being  monitored  also  shuts  down  immediately 
(for  example,  due  to  safety  reasons).  However,  in  many  situations  it  is  not  safe  or 
desirable  for  the  system  being  monitored  to  be  shut  down  immediately. 

In  situations  where  the  system  does  not  shut  down  immediately,  there  are  strong 
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reasons  why  one  may  desire  that  the  OVE  and  OVETV  algorithms  continue  to  track 
the  parameters  moving  in  the  parameter  space.  If  the  algorithm  is  being  used  in  an 
indirect  adaptive  control  setting,  the  parameter  set  is  critical  for  the  on-line  control 
design.  To  isolate  the  type  of  fault  that  has  occurred,  it  is  important  that  we  know 
where  in  the  parameter  space  the  parameter  set  has  moved  to  after  a  fault  has  been 
detected.  In  [48],  it  is  argued  that  the  need  for  estimation  of  uncertainty  in  the  pa¬ 
rameter  estimates  (which  is  integral  to  OVE  and  OVETV  algorithms)  is  particularly 
important  for  reconfigurable  control.  Fault  isolation  will  be  discussed  further  in  our 
concluding  remarks. 

In  [47],  a  recovery  strategy  is  developed  for  the  OBE  algorithm  of  [3].  In  [14],  a 
recovery  strategy  is  developed  for  a  modified  version  of  the  OBE  algorithm  [4].  When 
an  inconsistency  is  detected,  they  can  guarantee,  under  certain  assumptions,  that  this 
algorithm  will  asymptotically  contain  the  “true”  parameter. 

However,  we  would  argue  that  while  asymptotically  capturing  the  “true”  param¬ 
eter  is  a  nice  property,  it  does  not  satisfy  the  requirement  that  the  “true”  parameter 
always  be  contained  within  the  parameter  set.  We  suggest  instead,  that  after  a  fault 
is  detected,  the  ellipsoid,  Ek,  be  reset  to  a  “large  enough”  volume  ellipsoid  such  that 
the  “true”  parameter  is  contained  within  the  parameter  set. 

However,  initially  this  reset  ellipsoid  may  be  too  large  for  effective  robust  adaptive 
control,  isolation,  etc.  Consequently,  an  alternative  approach,  which  uses  projections 
to  effectively  utilize  the  information  gained  by  new  measurement,  Fk,  and  the  infor¬ 
mation  contained  in  the  ellipsoid,  Ek,  is  developed.  Finally,  an  integrated  approach  is 
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proposed  to  take  advantage  of  the  guaranteed  containment  properties  of  the  resetting 
algorithm  and  the  improved  volume  properties  of  the  projection  algorithm. 

5.4.1  Ellipsoid  Resetting  Algorithm 

Consider  the  example  investigated  in  Section  5.3.1.  This  is  similar  to  the  first-order 
example  considered  in  Section  3.2.  However,  at  A:  =  16  a  fault  occurs  causing  the 
parameter  to  jump  by  0.25,  i.e.,  Oie  =  flis  +  0.25.  A  simulation  was  run  and  the  fault 
was  detected  at  time  k  =  16;  the  intersection  of  Eiq  and  Fiq  was  empty  as  can  be 
seen  in  Figure  34.  At  this  point  the  OVETV  algorithm  halts.  As  discussed  above, 
often  we  would  like  to  continue  to  track  the  parameters  in  the  parameter  space. 

In  this  section,  we  will  reset  the  ellipsoid  to  a  size  which  is  guaranteed  to  capture 
the  true  parameters  after  the  fault  has  been  detected.  In  the  example  above,  the 
ellipsoid  Eiq  was  reset  to  have  Ae  =  27  and  =  0.  The  results  are  shown  in 
Figure  37.  The  algorithm  tracks  the  jump  in  the  parameter  very  quickly.  In  fact,  it 
appears  that  the  actual  parameter  is  always  contained  in  the  parameter  set. 

This  is  verified  by  examining  the  normalized  center  estimate  error  (see  Figure  38), 
which  is  found  to  always  be  less  than  one.  While  the  OVETV  algorithm  is  guaran¬ 
teed  to  contain  the  true  parameter  before  a  fault  occurs,  and  the  ellipsoid  resetting 
algorithm  guarantees  this  property  immediately  after  a  fault  is  detected,  there  is 
no  guarantee  that  the  true  parameter  will  always  be  contained  during  the  time  be¬ 
tween  when  the  fault  occurs  and  when  it  is  detected.  For  this  example,  the  fault  was 
detected  immediately;  this  will  not  always  be  the  case. 
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Figure  37:  Ellipsoid  resetting  algorithm  results  for  First-Order  example 
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Figure  38:  Normalized  center  estimate  error  for  First-Order  example 
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Furthermore,  the  ellipsoid  resetting  algorithm  guarantees  that  the  true  parame¬ 
ter  will  always  be  contained  immediately  after  a  fault  is  detected.  Under  nominal 
conditions,  once  the  parameter  is  captured,  it  is  guaranteed  to  remain  inside  the  pa¬ 
rameter  set  until  another  fault  occurs.  However,  if  the  system  does  not  satisfy  the 
nominal  conditions  after  a  fault  is  detected,  e.g.  due  to  continued  large  jumps  in  the 
parameter  vector,  the  true  parameter  may  not  remain  in  the  parameter  set.  Hope¬ 
fully,  however,  another  inconsistency  will  indicate  the  “continued”  fault  allowing  the 
true  parameter  to  be  recaptured  permanently.  This  is  an  important  point  and  will 
be  discussed  further  in  Section  5.5.1. 

One  of  our  concerns  about  the  ellipsoid  resetting  algorithm  was  the  volume  of  the 
ellipsoids  immediately  after  the  resetting  tahes  place.  Indeed,  as  seen  in  Figure  39,  the 
ellipsoid  volume  does  rise  sharply  immediately  after  the  ellipsoid  is  reset.  However, 
due  to  the  excitation  of  the  input  and  the  low  order  of  the  problem,  the  ellipsoid  vol¬ 
ume  decreases  rapidly.  In  the  next  section,  we  will  examine  an  alternative  algorithm 
which  should  result  in  smaller  ellipsoid  volumes,  but  lacks  the  guaranteed  recapture 
property  of  the  ellipsoid  resetting  algorithm. 

The  procedure  for  incorporating  the  resetting  algorithm  within  the  standard  OVETV 
structure  is  given  as  follows: 

1.  Choose  the  initial  ellipsoid,  Eo,  “large  enough”  such  that  the  initial  “true” 
parameter  vector,  6o,  is  in  Eq. 

2.  Find  Ek  such  that 


Ek  =  arg£  min{vol(Fl)  :  E  D  Gk}- 


(5.32) 
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Figure  39:  Ellipsoid  volume  for  First-Order  example 


3.  If  Efc  n  Ffe  =  0,  then  signal  a  fault  and  reset  Ek, 

Ek  —  Ereseti  (5.33) 

where  Ereset  is  “large  enough”  such  that  “true”  parameter  vector,  9k,  is  in  Ek- 

4.  Find  Ek  such  that 

Ek  =  arg£  min{vol(F)  :  E  D  Ek(^  Ffc}.  (5.34) 

5.  Repeat  steps  two  through  four  for  each  new  measurement. 

Note,  that  in  step  3,  the  question  of  how  one  chooses  Ereset  “large  enough”  such  that 
9k  ^  Ek  is  an  important  one  and  will  be  discussed  in  detail  in  Section  5.5.2. 
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5.4.2  Ellipsoid  Projection  Algorithm 

Consider  the  situation  in  Figure  40  where,  at  time  k,  the  “true”  parameter  6k  has 
jumped.  The  intersection  of  Fk  and  Ek  is  empty,  and  a  fault  has  been  detected  by  the 
OVETV  algorithm.  In  this  section,  we  want  to  develop  an  algorithm  to  find  a  new 
ellipsoid,  Ek,  which  utilizes  the  information  gained  by  Fk  and  the  “best”  information 
possessed  by  Ek. 

We  are  guaranteed  that  the  true  parameter  must  be  contained  within  Fk.  The 
region,  Fk,  provides  information  in  the  direction  parallel  to  the  regression  vector,  <f)k. 
However,  it  provides  no  information  in  the  directions  orthogonal  to  (j)k.  Since  we 
have  no  additional  information,  we  will  make  the  assumption  that  in  the  directions 
orthogonal  to  (f>k,  that  the  ellipsoid,  Ek,  bounds  the  “true”  parameter.  While  this 
assumption  holds  for  the  example  in  Figure  40,  there  is  no  guarantee  that  it  will  hold 
in  general.  We  will  use  the  results  of  Sections  2.4.2  (optimal  time  update  equations) 
and  4.3  (dependent  parameter  variations)  to  solve  for  the  minimum  volume  ellipsoid 
which  bounds  the  sum  of  (i)  the  degenerate  ellipsoid  found  by  projecting  Ek  onto  the 
plane  orthogonal  to  ^k  and  (ii)  the  degenerate  ellipsoid  found  by  projecting  Fk  onto 
the  line  of  (f>k.  In  Figure  40,  the  resulting  ellipsoid  is  given  by  Ek. 

Mathematically,  we  want  to  find  Ek  such  that 

Ek  =  arg£;min{vol(E)  :  E  D  V{Ek,^i)  +V{Fk,(f>k)},  (5.35) 

where  V{X,  Y)  is  the  projection  of  X  onto  Y  and  Y-^  is  the  orthogonal  complement 
of  Y.  To  solve  this  problem,  we  begin  by  finding  the  singular  value  decomposition  of 
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Figure  40:  Projection  Algorithm 


<f>k  = 


(5.36) 


where  U  G  is  orthogonal,  S  G  is  equal  to  S  =  cr  0  •  •  •  0  ]  ,  and  V  G  3^ 
is  equal  to  one.  It  is  easy  to  show  that  a  must  be  equal  to  ||^fc||2-  Recall  that  Fk  and 
Ek  are  defined  as 


Fk  =  {0fe€r  (5.37) 

Ek  =  {Ok  :  {Ok  -  hfPk\0k  -  k)  <l\0ke  r },  (5.38) 


where  h  €  is  the  center  of  the  ellipsoid  and  Pk  €  is  a  symmetric,  positive 


definite  matrix. 
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To  simplify  our  problem,  we  define  the  affine  transformation,  Tu,  and  its  inverse 

by 

h  =  Tu{0k)  =  U^{0k-h)  (5.39) 

0k  =  T{j\0k)  =  U0k  +  h^  (5.40) 

This  transformation  is  similar  in  form  to  the  transformation  defined  in  [29]  for  the 
development  of  the  MOVE  algorithm.  The  key  difference  lies  in  the  fact  that  U  is 
orthogonal  and  therefore  directions  and  distances  are  preserved. 

Substituting  (5.40)  into  (5.38),  we  see  that  the  ellipsoid,  Ek,  is  given  by 

Ek  =  {0k  :  0lPT"h  <  1;4  G  r}  (5.41) 

in  the  transformed  coordinate  systems  where  Pk  =  U^PkU.  The  transformation 
translates  Ek  to  the  origin  and  then  applies  a  rotation.  Substituting  (5.40)  into 
(5.37),  we  see  that  the  region  between  the  two  hyperplanes,  Fk,  is  given  by 

h  =  [h  e  :  7^  <  j/fc  -  -  0l(l>k  <  Ik)  (5.42) 

in  the  transformed  coordinate  systems  where  =  U'^ <f>k-  Using  (5.36),  ^k  can  be 
rewritten  as 

^k  =  =  U^UEV'^  =  S.  (5.43) 

Substituting  (5.43)  into  (5.42)  and  rearranging  terms  results  in 

Fk  =  {0k  €  :  Ofc  <  0ki  <  ajt},  (5.44) 


where  ak  and  ock  are  now  defined  as 
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Oik 


Vk  -  n^k-i 


±k 


Ukh 


(5.46) 


and  9k  =  ^fc2  •  •  •  9kr  ]  •  The  variables  ^  and  ak  represent  the  location  of 

the  parallel  hyperplanes  which  are  orthogonal  to 
Let  the  matrices,  Pk  and  U,  be  partitioned  as 

C/f 


Pk  = 


Ai 

A2' 

[  Ai 

P22 

U^  = 


uT 


(5.47) 


where  Ai  G  3?,  P12  G  Ai  G  A2  G  and 

U2  G  The  projection  of  Ek  onto  the  subspace  orthogonal  to  ^k  is  given  by 

the  degenerate  ellipsoid 


V{Ek,  ^i)  =  {h  G  r  :  rjj9k  <  ^/riJQiVr]  Vr  G  ^ (5.48) 


where 


Qi  = 


■  0 

0 

0 

to 

(5.49) 


The  projection  of  Fk  onto  the  axis  parallel  to  is  given  by  the  degenerate  ellipsoid 
V{Fk,h)  =  {Ok  G  r  :  vJOk  <  rjjh  +  G  r},  (5.50) 

where 


Q2  — 

and  a  =  ^  and  /?  = 


0  ■ 

0 

0 

a 


(5.51) 


Next,  we  want  to  solve  (5.35)  in  the  transformed  coordinate  systems,  i.e., 


Ek  =  argjE  min{vol(£)  ;  E  D  V{Ek,  ^i)  +  V{Fk,  ^k)]- 


(5.52) 
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This  is  the  problem  we  solved  in  Section  2.4.2  for  the  case  where  both  ellipsoids  were 
nondegenerate  and  in  Section  4.3  for  the  case  where  one  ellipsoid  was  nondegenerate. 
However,  in  (5.52)  both  ellipsoids  are  degenerate.  To  apply  the  results  of  Section  2.4.2, 
we  will  let 


e 

0 

0 

P22 

(5.53) 


where  e  is  a  scalar  which  we  will  specify  later,  and  factor  Q2  as 


Q2  =  BRB^, 


(5.54) 


where  B^  = 


10 


and  R 


From  Section  2.4.2  we  know  that  the  solution  of  (5.52)  can  be  found  by  solving 
the  following  generalized  eigenvalue  problem 


BRB^Xj  —  \jQiXj. 


(5.55) 


for  each  Aj,  j  G  [l,r].  However,  from  (4.34)  in  Section  4.3,  we  know  that  the  solution 
of  (5.52)  can  also  be  found  by  the  solving  the  following  reduced  order  problem 

B^Q^^Bx  =  \R~^x  (5.56) 

for  A.  This  is  equivalent  to  solving 

det(Ai2-^  -  B^Q^^B)  =  0  (5.57) 

for  A.  The  solution  to  (5.57)  is  given  by 
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For  the  case  where  s  =  1,  the  unique  value  of  >  0  in  (4.35)  is  given  by 


P  = 


-A(r  —  1)  +  ^A2(r  —  1)^  +  4rA 


2r 


(5.59) 


After  some  algebra,  (5.59)  can  be  rewritten  as 

2 


P  = 


(5.60) 


Substituting  (5.58)  into  (5.60)  and  taking  the  limit  as  e  approaches  0,  we  find  that 


P  = 


r  —  1 


(6.61) 


Therefore,  from  (4.36)  and  (4.37),  we  find  that  the  solution  to  (5.52)  is  given  by 

E,  =  {4  :  (h  -  -  4)  <  1;  4  e  r },  (5.62) 

where 


Pk  =  {p  +  1)^2  +  (P"'  +  = 


0 

0 

J 

(5.63) 


Substituting  (5.39)  into  (5.62)  gives  us  the  solution  to  (5.35)  in  the  original  coor¬ 
dinate  system  as 


Ek  =  {Ok :  {0k  -  0kyPk\0k  -  0k)  <  1;  4  e  r }, 


(5.64) 


where 


Pfc  =  ry'^UiU'^  + —-U2P22U2 


r  —  1 


4  —  4  +  olU\. 


(5.65) 

(5.66) 
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By  substituting  for  o;  and  j3  and  using  the  fact  that  4>k  —  crUi,  (5.65)  and  (5.66)  can 
be  rewritten  as 

0k  =  0k +  - - ^Vt - —<l>k-  (5.68) 


Therefore,  the  ellipsoid  projection  algorithm  for  recovering  from  a  fault  is  given 
as  follows: 

1.  Find  the  singular  value  decomposition  of  <f>k  as 


cf>k  =  UEV^. 

(5.69) 

2.  Calculate 

Pk  =  U'^PkU. 

(5.70) 

3.  Update  0k  and  Pk 

Vk  <i>l0k 

6k  =  ■■  'A — —<i>k 

<Pk  9k 

(5.71) 

5  _  ,  ,T.  TT  p  JjT 

(5.72) 

As  discussed  earlier,  one  iteration  of  the  projection  algorithm  is 

not  guaranteed 

to  capture  the  “true”  parameter.  Consequently,  even  after  the  projection  algorithm 
is  applied  and  the  OVETV  algorithm  resumes,  more  inconsistencies  could  arise.  If 
this  happens,  the  projection  algorithm  would  then  be  applied  again. 

A  procedure  for  incorporating  the  projection  algorithm  within  the  standard  OVETV 


structure  is  given  as  follows: 
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1.  Choose  the  initial  ellipsoid,  Eq,  “large  enough”  such  that  the  initial  “true” 
parameter  vector,  Oq,  is  in  Eq. 

2.  Find  Ek  such  that 

.Efc  =  arg£;min{vol(£;)  :  D  Gfc}.  (5.73) 


3.  Find  Ek  such  that 

EkC\Fk  if  n  Ffc  7^  0  ^ 

V{Ek,(f>k)  +'P{Ek,(l)k)  otherwise 

(5.74) 

4.  Repeat  steps  two  and  three  for  each  new  measurement. 


Ek  =  arg£;min< 


vol(F?)  :  E  D 


This  procedure  is  often  able  to  recapture  the  “true”  parameter  very  quickly.  Con¬ 
sider  the  first-order  example  we  investigated  in  Section  5.4.1  where  a  fault  occurs  at 
A:  =  16  causing  the  ak  parameter  to  jump  by  0.25,  i.e.,  aie  =  ai5  -f  0.25.  This  simula¬ 
tion  was  repeated  using  the  projection  algorithm,  rather  than  the  resetting  algorithm, 
as  a  recovery  strategy.  The  results  are  shown  in  Figure  41. 

For  this  example,  the  results  are  very  similar  to  those  obtained  using  the  resetting 
algorithm.  The  key  difference  can  be  seen  in  the  ellipsoid  volumes  of  Figure  39.  The 
large  increase  in  ellipsoid  volume,  which  occurs  immediately  after  the  fault  is  detected 
in  the  resetting  algorithm  results,  is  not  seen  in  the  projection  algorithm  results.  As 
can  be  seen  in  Figure  38,  the  true  parameter  is  always  contained  in  the  parameter 
set  for  both  simulations.  This  is  particularly  surprising  for  the  projection  algorithm, 
where  repeated  applications  of  the  algorithm  are  often  necessary  to  recapture  the 
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Figure  41:  Ellipsoid  projection  algorithm  results  for  First-Order  example 


true  parameter.  In  the  next  section,  we  will  suggest  a  set  of  rules  which  serve  to 
combine  the  ellipsoid  resetting  algorithm  and  the  ellipsoid  projection  algorithm  in  an 
integrated  approach. 


5.4.3  Integrated  Approach 

The  ellipsoid  resetting  algorithm  is  guaranteed  to  capture  the  true  parameter  imme¬ 
diately  after  a  fault  is  detected.  This  is  important  for  robust  control  applications  and 
detection  of  any  future  faults.  However,  immediately  after  being  reset,  the  resulting 
ellipsoid  may  be  too  large  for  effective  control  applications.  On  the  other  hand,  the 
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projection  algorithm  does  not  suffer  from  the  large  volume  increases  encountered  by 
the  resetting  strategy.  While  often  capturing  the  true  parameter  soon  after  a  fault 
is  detected,  the  projection  algorithm  may  have  to  “recover”  from  several  inconsisten¬ 
cies  before  doing  so.  Furthermore,  it  lacks  the  guaranteed  properties  of  the  resetting 
algorithm. 

In  this  section,  we  propose  an  algorithm  whereby  once  a  fault  is  detected  by  an 
inconsistency,  the  resetting  algorithm  and  projection  algorithm  are  run  in  parallel.  A 
set  of  rules  are  used  to  decide  which  ellipsoid  to  utilize  for  robust  control  applications. 
Although  heuristic,  these  rules  are  designed  to  capture  the  strengths  of  each  approach. 
One  possible  set  of  rules  is  given  below,  where  TU  are  the  time  update  equations  of 
Section  2.4.2,  MU  are  the  measurement  update  equations  of  Section  2.3,  and  PU  are 
the  projection  update  equations  of  Section  5.4.2.  In  this  notation,  “MU(jFa;,  Fk)",  for 
example,  means  to  apply  the  measurement  update  equations  to  the  ellipsoid  Ek  and 
region  between  two  hyperplanes  Fk.  Once  a  fault  has  been  detected,  El  and  El  are 
the  ellipsoids  resulting  from  the  resetting  algorithm,  and  E^  and  El  are  the  ellipsoids 
resulting  from  the  projection  algorithm. 

The  steps  of  this  integrated  approach  are  given  below. 

1.  Choose  the  initial  ellipsoid,  Eq 

2.  k  =  kAl 

3.  Ek —  A:\}{Ek-i.,Rk-i) 


4.  If  EkV\Fk^  0,  then 
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.  Ek  =  M\]{Ek,Fk) 

•  Go  to  step  2 

5.  A  fault  has  been  detected:  E^  =  PU(^fc,  4k,  Fk)  El  =  M\J{Ereset,  Fk) 

6.  If  (a)  EkC\Ek  =  0,  (b)  vol(£Jfc)  <  vol(£^^),  or  (c)  vol(£^^)  <  7,  then 

•  Ek  =  El 

•  Go  to  step  2 

7.  Ek  =  El 

8.  k  =  k  +  1 

9.  El  =  T\]{El_,,Rk.^)  El  =  T\]{EU,Rk-i) 

10.  If  EinFk  =  0,  then 

•  Ek  =  El 

•  Go  to  step  5 

11.  El  =  ME{ElFk) 

12.  If  ElnFk  =  0,  then 

.  El  =  P\]{El,4k,Fk) 

•  Go  to  step  6 

13.  El  =  MV{P„Fi,) 
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14.  Go  to  step  6 

Of  particular  importance  in  this  integrated  approach  are  steps  6  and  10.  Once  a 
fault  has  been  detected,  they  are  used  to  decide  whether  the  ellipsoids  resulting  from 
the  projection  strategy  or  the  ellipsoids  resulting  from  the  resetting  strategy  are  to  be 
used  for  control  purposes.  In  step  6,  rule  (a)  is  the  most  important.  If  E^nEl  —  0,  El 
does  not  contain  the  true  parameter;  therefore,  we  should  use  El  for  robust  control. 
In  rule  (b),  we  choose  El  if  it  has  the  smallest  volume  because  it  is  guaranteed  to 
contain  the  true  parameter,  while  the  projection  strategy  is  not.  In  rule  (c),  we  also 
choose  El  if  the  volume  is  “small  enough”  that  a  robust  controller  design  is  possible, 
where  7  is  a  prespecified  constant.  Clearly,  other  rules  could  be  used  instead  of  (b) 
and  (c)  depending  on  the  control  design,  performance  requirements,  etc.. 

Step  10  is  used  to  determine  if  an  additional  fault  has  occurred.  Notice  that  the 
determination  of  a  fault  depends  on  El  and  not  on  El  because  of  the  guaranteed 
containment  property  of  the  resetting  algorithm. 

Probably  the  biggest  disadvantage  of  the  integrated  approach  is  the  added  com¬ 
plexity.  After  recovery,  and  until  El  has  been  selected,  two  versions  of  the  OVETV 
algorithm  are  required  to  run  in  parallel. 


5.5  Example 

In  this  section,  we  will  consider  the  linear  time- varying  circuit  which  was  examined  in 
Section  3.3.  Recall  that  we  investigated  a  series  R-L-C  circuit  where  the  components 
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were  allowed  to  vary  with  time.  To  obtain  the  true  parameters,  the  continuous  time 
state-space  equations  were  discretized  and  then  transformed  into  a  difference  equation 
using  Theorem  5.  Quantization  errors  in  the  D/A  and  A/D  converters  were  also 
incorporated  into  the  ARX  structure. 

In  the  simulations  of  Section  3.3,  the  inductance  I  and  capacitance  c  were  held 
constant,  while  the  resistance  r  was  allowed  to  vary  slowly  with  time.  In  this  section, 
we  will  consider  two  scenarios.  The  values  of  I  and  c  will  remain  constant  for  both 
scenarios.  In  the  first  scenario,  the  nominal  model  will  assume  that  r  varies  slowly 
with  time.  The  consistency  check,  which  is  integral  to  the  OVETV  algorithm,  will 
be  used  to  detect  when  an  abrupt  fault  occurs.  The  ellipsoid  resetting  and  ellipsoid 
projection  algorithms  will  be  applied  to  recover  tracking. 

In  the  second  scenario,  the  system  is  nominally  time-invariant.  Our  goal  is  to 
detect  incipient  faults  and  track  the  parameters  after  a  fault  is  detected.  Two  different 
methods  are  used  to  detect  a  fault.  In  the  first  method,  the  OVETV  algorithm  is 
applied  to  the  system  with  the  bounds  on  the  parameter  disturbance  vector  set  to 
track  the  parameters  during  an  incipient  fault.  The  ellipsoid  intersection  test  detects 
a  fault  when  the  ellipsoid  output  of  the  OVETV  algorithm  does  not  intersect  with 
the  ellipsoid  which  was  found  using  the  OVE  algorithm  under  nominal  conditions. 
In  the  second  method,  the  OVE  algorithm  is  applied  to  the  system  and  the  internal 
consistency  check  is  used  to  detect  a  fault.  The  resetting  and  projection  schemes  will 
be  applied  to  recover  tracking. 
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5.5.1  Abrupt  Fault 

In  this  section,  we  will  use  the  consistency  check  of  the  OVETV  algorithm  to  detect 
an  abrupt  fault  in  the  linear  LTV  circuit  of  Section  3.3.  Nominally,  we  expect  the 
circuit  to  behave  just  as  it  did  in  Section  3.3.  Consequently,  the  OVETV  algorithm 
is  initialized  the  same  with  j).  =  — 7^,  =  0.005,  Rk  —  (ll  where  (k  =  0.001,  Po  =  107, 
and  00  =  0.  The  input,  Umk,  is  a  uniformly  distributed  random  variable  between 
[-5.0, 5.0]  volts,  and  the  sampling  time  T  is  0.1  seconds. 

In  Section  3.3,  the  actual  circuit  components  were:  r{t)  =  (1.0  +  0.5cos(t/10.0))0, 
l{t)  =  l.OH,  and  c{t)  =  l.OF.  In  this  section,  the  components  are  identical  to  those  of 
Section  3.3  from  t  =  0.0  to  t  =  10.0.  However,  at  t  =  10.0,  the  resistance  r{t)  jumps 
sharply,  i.e.,  r{t)  =  (1.0  +  0.5  cos(t/10.0)  +  O.busit  -  10))O,  where  U5(-)  is  the  unit 
step  function. 

A  simulation  was  run  with  the  OVETV  consistency  check  utilized  for  detection 
and  the  ellipsoid  projection  scheme  utilized  for  recovery.  The  consistency  check  easily 
detected  the  fault  and  signalled  immediately  at  t  =  10.0  seconds. 

Consequently,  the  projection  scheme  was  applied  to  recover  tracking.  Inconsisten¬ 
cies  also  occurred  at  t  =  10.2  and  t  =  10.3  seconds  requiring  a  reapplication  of  the 
projection  scheme.  This  reapplication  was  expected  since  the  projection  algorithm  is 
not  guaranteed  to  recapture  the  parameter.  Finally,  as  seen  in  Figure  42,  the  normal¬ 
ized  center  estimate  error  goes  below  1.0  at  t  =  10.6  seconds.  This  tells  us  that  the 
true  parameter  has  been  recaptured  by  the  parameter  set.  Under  nominal  conditions, 
once  the  true  parameter  has  been  recaptured,  it  will  stay  there  until  another  fault 
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Seconds 

Figure  42:  Normalized  center  estimate  error  during  abrupt  fault 

occurs.  Figure  43  shows  the  tracking  of  the  parameter  bi  before  and  after  the  fault 
occurs. 

The  simulation  was  run  again;  this  time  using  the  ellipsoid  resetting  scheme  for 
recovery.  The  resetting  ellipsoid,  Preset was  set  the  same  as  Eq,  i.e,  Preset  =  10/ 
and  9reset  =  0.  As  before,  the  fault  was  detected  immediately  at  t  —  10.0  seconds. 
As  expected,  the  ellipsoid  resetting  strategy  resets  Ek  to  Ereset  before  applying  the 
measurement  update  equations.  As  can  be  seen  in  Figure  42,  the  normalized  center 
estimate  error  for  the  ellipsoid  resetting  algorithm  is  clearly  less  than  one  att  =  10.0 
seconds. 

At  this  point,  note  that  the  OVETV  signals  another  inconsistency  at  t  =  10.4 
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Figure  43:  Parameter  bi  with  Ellipsoid  Projection  algorithm 


seconds.  Intuitively,  this  may  be  unexpected  since  the  ellipsoid  resetting  algorithm 
was  designed  to  recapture  the  parameter  immediately  after  a  jump  in  the  parameter 
vector  was  detected.  While  the  algorithm  succeeds  in  doing  this  at  i  =  10.0  and 
t  =  10.4  seconds,  as  can  be  seen  in  Figure  42,  the  normalized  center  estimate  error  is 
greater  than  one  between  these  times. 

To  understand  these  results,  consider  the  change  in  the  true  parameter  vector  at 
each  time  step,  Wk.  Under  nominal  conditions,  Wk  is  bounded  by  (2.4).  Recall  from 
above  that  Rk  was  specified  to  be  equal  to  where  (fc  =  0.001.  When  Rk  =  (k^^ 
the  bound  in  (2.4)  can  be  written  as 


lltyfcih  <  Cfc- 


(5.75) 
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Figure  44:  Change  in  the  true  parameter  vector 


In  Figure  44,  ||t0fc||2  and  (k  are  plotted.  As  expected,  at  t  =  10.0,  ||r«A:||2  >  Cfc- 
Clearly,  a  fault  has  occurred.  However,  at  t  =  10.1  and  t  =  10.2,  l|tufc||2  is  again 
greater  that  (k-  This  is  puzzling  at  first,  since  the  only  jump  in  the  continuous  time 
state  space  equations  occurred  at  t  =  10.0.  However  by  examining  equation  (3.32), 
we  see  that  true  parameter  vector,  Ok  is  a  function  of  ^2  at  time  k,  k  —  1,  and  k  —  2, 
where  d2{k)  =  r{kT). 

Because  the  OVETV  algorithm  is  highly  dependent  on  the  input  of  the  system, 
we  ran  this  simulation  again  using  a  saw-tooth  wave  as  an  input.  The  results  were 
similar,  but  the  inconsistencies  were  detected  at  t  =  10.0  and  t  =  10.2. 

This  problem  poses  an  interesting  challenge  for  a  recovery  strategy,  because  even 
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Figure  45:  Parameter  bi  with  Ellipsoid  Resetting  algorithm 


though  the  jump  in  the  continuous  time  parameter  occurred  during  just  one  time 
step,  it  caused  the  nominal  system  bounds  to  be  violated  for  r  time  steps.  Perhaps 
the  integrated  approach,  which  was  discussed  in  Section  5.4.3,  could  be  modified  to 
account  for  this  behavior. 

Lest  we  forget,  this  algorithm  actually  performed  quite  well.  It  was  able  to  detect 
the  fault  immediately,  recapturing  the  true  parameter  one  time  step  after  the  fault 
ended.  Figure  43  shows  the  tracking  of  the  parameter  bi  before  and  after  the  fault 
occurs.  As  would  be  expected,  the  ellipsoid  projection  algorithm  had  a  smaller  volume 
around  the  fault  area.  See  Figure  46. 
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Figure  46:  Ellipsoid  volume  during  abrupt  fault 


5.5.2  Incipient  Fault 

In  the  example  of  the  previous  section,  the  nominal  plant  was  slowly  time-varying, 
but  an  abrupt  fault  occurred  in  the  actual  plant.  In  this  example,  the  nominal  plant 
is  time-invariant  and  the  fault  is  incipient  (slowly- varying).  Two  methods  will  be 
used  to  detect  the  fault. 

The  first  method  uses  the  ellipsoid  intersection  test  of  section  5.3.2  to  detect 
a  fault.  First,  we  apply  the  OVE  algorithm  to  the  system  when  it  is  operating 
under  nominal  conditions.  Under  nominal  conditions,  the  circuit  components  are 
r(t)  =  1.5fl,  l(t)  =  l.OH,  and  c(t)  =  l.OF.  The  OVE  algorithm  is  initialized  with 

=  —7^  =  0.005,  Po  =  107,  and  0o  =  0.  After  20  seconds,  the  bounding  ellipsoid  is 
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given  by  Enom,  where 


1.7476e  -  03 
-1.5366e  -  03 
2.6452e  -  06 
2.4247e  -  04 


— 1.5366e  —  03 
1.3697e  -  03 
-3.5244e  -  06 
-2.1448e  -  04 


2.6452e  -  06 
-3.5244e  -  06 
1.9640e  —  06 
5.4088e  -  07 


2.4247e  -  04  ' 
-2.1448e  -  04 
5.4088e  -  07 
3.5318e  -  05 


(5.76) 


Onrtm  — 


-1.8537e  +  00  8.6256e  -  01  1.3916e-01  -1.3935e  -  01  .  (5.77) 


When  the  system  is  undergoing  an  incipient  fault,  it  is  still  known  to  satisfy  the 
bounds  which  were  specified  in  Section  3.3.  Consequently,  to  monitor  the  circuit, 
the  OVETV  algorithm  is  initialized  the  same  as  in  Section  3.3.  The  system  was 
simulated  with  the  actual  circuit  components  given  by  r{t)  =  (1.0  +  0. 5  cos(t/10.0))fl, 
l{t)  =  l.OH,  and  c{t)  =  l.OF. 

The  ellipsoid  intersection  test  is  used  to  indicate  when  Ek,  the  output  of  the 
OVETV  algorithm,  and  Enom,  the  nominal  ellipsoid,  fail  to  intersect.  At  t  =  3.6 
seconds,  the  ellipsoids  fail  to  intersect  and  a  fault  is  indicated.  The  performance  of 
this  algorithm  clearly  depends  on  how  tightly  Ek  bounds  the  actual  parameter  and 
the  size  of  the  nominal  ellipsoid,  Enom-  For  comparison  purposes,  it  was  found  that 
the  true  parameter  actually  escaped  the  nominal  parameter  set  at  t  =  2.0.  As  we  saw 
in  Section  3.3,  the  OVETV  algorithm  has  no  problem  tracking  this  system. 

The  second  method  we  will  apply  uses  the  consistency  check  which  is  integral 
to  the  OVE  algorithm  to  indicate  a  fault.  The  OVE  algorithm  was  initialized  to 
bound  the  nominal  system  with  7^  =  —7^  =  0.005,  Po  =  10/,  and  Oq  =  0.  The 
OVE  algorithm  was  then  applied  to  the  measured  data  from  the  actual  system.  At 
t  =  3.9  seconds,  an  inconsistency  was  indicated,  signaling  a  fault.  This  was  slightly 
slower  than  the  results  achieved  using  the  intersection  test.  However,  this  algorithm 
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Figure  47:  Normalized  center  estimate  error  during  incipient  fault 


is  significantly  cheaper  in  terms  of  the  number  of  computations  required  at  each  time 
step.  It  requires  neither  the  time  update  computations  of  the  OVETV  algorithm  nor 
the  computations  required  for  the  ellipsoid  intersection  test. 

The  ellipsoid  intersection  strategy  has  the  advantage  of  built-in  tracking  after  an 
incipient  fault  is  detected.  This  raises  the  question  of  how  well  the  OVE  algorithm 
can  track  the  system  after  a  fault  is  detected  using  the  projection  and  resetting 
algorithms.  With  the  resetting  ellipsoid  initialized  to  Oreset  =  ^  and  Preset  =  10007, 
the  normalized  center  estimate  errors  for  the  OVE  algorithm  with  ellipsoid  resetting 
and  ellipsoid  projection,  and  for  the  OVETV  algorithm,  are  given  in  Figure  47. 

Clearly,  the  OVETV  algorithm  has  the  best  tracking  characteristics  in  terms  of 
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Table  6:  Center  estimate  performance  measures 


OVETV 

Ellipsoid 

Resetting 

Ellipsoid 

Projection 

£2 

0.0800 

0.1019 

0.1358 

0.1780 

0.0589 

0.0740 

normalized  center  estimate  error,  as  it  always  remains  below  one.  The  ellipsoid  pro¬ 
jection  normalized  center  estimate  error  resets  below  one  every  time  an  inconsistency 
occurs,  but  drifts  above  one  as  the  parameters  continue  to  vary  with  time.  The  results 
for  the  ellipsoid  projection  algorithm  are  clearly  unacceptable  since  they  spend  much 
of  the  time  above  one. 

While  our  primary  interest  is  in  the  parameter  sets,  it  is  interesting  to  note  that  the 
center  estimate  of  the  ellipsoid  projection  algorithm  actually  tracks  the  true  parameter 
better  than  either  of  the  other  approaches.  Based  on  the  performance  measures 
defined  in  Section  3.3,  the  projection  algorithm  results  have  the  lowest  cost  as  can  be 
seen  in  Table  6. 

Having  gained  some  experience  with  the  ellipsoid  projection  algorithm,  this  seems 
to  be  an  appropriate  place  to  discuss  how  Ereset  can  be  chosen  “large  enough”.  The 
answer  is  definitely  problem  dependent.  If  one  has  a  relatively  significant  amount  of 
a  priori  information  about  the  system,  i.e.,  how  it  behaves  under  nominal  and  fault 
modes,  then  it  may  be  possible  to  specify  Ereset  for  the  particular  case. 

If,  on  the  other  hand,  relatively  little  is  known  about  the  system,  at  least  two 
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choices  exist.  In  the  abrupt  fault  example,  we  chose  an  Ere$et  of  the  form  $reset  = 
0,  Preset  =  /5/  where  ^  is  a  large  number.  This  approach  basically  discounts  any 
measurements  received  prior  to  that  point. 

In  the  incipient  fault  example,  we  chose  an  Ereset  of  the  form,  Preset  =  &  and 
Preset  =  ^I  where  is  a  large  number.  This  gives  some  credibility  to  the  center 
estimate  of  the  previous  ellipsoid,  but  greatly  increases  the  uncertainty.  Observations 
made  using  this  form  for  Ereset  led  to  the  following  Theorem. 

Theorem  9  Let 

E{6r,Pr)  -  arg£;min{TO/(E)  :  E  D  E{9reset,  Preset)  H  Fk} ,  (5.78) 

and 

E{9p,  Pp)  =  avgE  min{vol{E)  :  E  D  V{E{h,  Pk),  ^i)  +  ViFk,  (i>k)},  (5.79) 

where  E{9,P)  is  defined  to  be 

E{9,P)  =  {9k  :  {9k  -  9)'^p-\9k  -  9)  <  1;4  e  r}.  (5.80) 

Given  that  E{9k,Pk)  C\  Fk  =  0,  if  9 reset  =  9k  and  Preset  —  fil,  then 

lim  9r  =  9p.  (5.81) 

P-*oo 

Proof:  The  proof  is  straight  forward.  From  (5.71), 


5.6.  Summary 


155 


It  is  easy  to  show  that  in  the  limit,  as  ^  approaches  infinity,  Or  equals  (5.82).  This 
is  done  by  substituting  Oreset  =  4  and  Preset  =1^1  into  the  measurement  update 
equations  (2.14)-(2.25)  of  Section  2.3  and  taking  the  limit  as  (3  approaches  infinity. 
□ 

This  theorem  is  significant  because  it  states  that  if  you  choose  the  parameters  of 
Preset  as  Oreset  =  Ok  and  Preset  =  where  j3  is  some  large  number,  then  the  center 
estimate  produced  by  the  resetting  algorithm  is  approximately  equal  to  the  center 
estimate  produced  by  the  projection  algorithm.  The  orientation  matrices,  however, 
are  quite  different. 


5.6  Summary 

In  this  chapter,  we  demonstrated  how  the  OVE  and  OVETV  algorithms  could  be 
used  for  EDI.  We  began  the  chapter  by  discussing  some  of  the  issues  involved  in  the 
implementation  of  EDI  algorithms. 

Next,  we  discussed  two  different  methods  for  detection  of  failures.  The  first 
method  used  the  consistency  check  which  is  integral  to  the  OVE  and  OVETV  al¬ 
gorithms.  Because  the  check  is  an  integral  part  of  the  OVE  and  OVETV  algorithms, 
this  detection  strategy  was  easy  to  implement.  It  was  particularly  adept  at  detecting 
abrupt  failures  in  nominally  slowly- varying  systems. 

We  also  utilized  an  ellipsoid  intersection  test  to  detect  when  an  ellipsoid,  produced 
under  nominal  conditions,  fails  to  intersect  with  the  current  ellipsoid  of  the  OVETV 
algorithm.  The  advantage  of  this  algorithm  lies  in  the  fact  that  it  is  able  to  track 
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incipient  faults  without  resorting  to  recovery  strategies.  However,  it  requires  addi¬ 
tional  computations  and  is  not  easily  applicable  to  nominally  time-varying  systems. 
An  alternative  development  of  the  ellipsoid  intersection  test  equations  was  realized 
by  framing  the  problem  as  a  constrained  optimization  problem. 

We  should  mention  that  the  ellipsoid  intersection  test  could  also  be  used  for  fault 
isolation.  Assume  that  “signature”  ellipsoids  exist  which  represent  the  system  un¬ 
der  particular  faults.  Once  a  fault  has  been  detected  and  recaptured,  the  ellipsoid 
intersection  test  could  be  used  to  check  which,  if  any,  signature  ellipsoids  intersect 
with  the  current  ellipsoid.  If  the  current  ellipsoid  does  not  intersect  with  a  particular 
signature  ellipsoid,  then  it  is  clear  that  the  associated  fault  did  not  occur.  Clearly, 
if  a  signature  ellipsoid  intersects  with  the  current  ellipsoid,  the  associated  fault  may 
have  occurred.  For  certain  systems,  it  may  be  possible  to  use  a  modified  version  of 
the  OVE-ISP  algorithm  in  Section  4.5  to  help  isolate  which  fault  did  occur. 

When  a  fault  is  detected  by  an  inconsistency  in  the  OVE  and  OVETV  algorithms, 
the  algorithms  halt.  For  applications  such  as  fault  isolation  and  reconfigurable  control, 
it  would  be  desirable  to  continue  tracking  the  parameters  after  the  fault  is  detected. 
Two  methods  for  recovering  the  true  parameter  were  proposed. 

The  ellipsoid  resetting  algorithm  resets  the  current  ellipsoid  to  one  large  enough 
to  capture  the  true  parameters.  If  the  system  satisfies  the  nominal  conditions  after 
the  ellipsoid  is  reset,  the  algorithm  is  guaranteed  to  contain  the  true  parameters  until 
another  fault  occurs.  Several  different  strategies  were  considered  for  choosing  the 
resetting  ellipsoid,  Ereset-  H  was  shown  that  for  a  particular  form  of  Ereseti  that  as 
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the  size  of  Ereset  is  increased,  the  center  estimate  approaches  the  center  estimate  of 
our  other  recovery  strategy,  ellipsoid  projection. 

To  guarantee  that  the  true  parameter  is  contained  in  the  parameter  set  after  a 
fault  is  detected,  the  ellipsoid  resetting  algorithm  often  results  with  large  ellipsoid 
volumes  immediately  after  a  fault  is  detected.  Consequently,  a  method  with  better 
volume  properties  for  recovering  the  true  parameter  is  also  proposed.  It  combines 
the  information  contained  in  the  new  measurements  and  assumes  correct  information 
from  the  old  ellipsoid  whenever  ambiguity  arises.  Because  this  is  not  always  the 
case,  repeated  applications  of  the  projection  algorithm  may  be  required  to  recapture 
the  true  parameter.  To  take  advantage  of  the  guaranteed  containment  properties  of 
the  ellipsoid  resetting  algorithm  and  the  smaller  volumes  of  the  ellipsoid  projection 
algorithm,  an  integrated  approach  was  also  proposed. 

Finally,  the  strategies  of  this  chapter  were  demonstrated  on  the  LTV  circuit  of 
Section  3.3.  The  algorithms  were  shown  to  successfully  detect  and  track  both  abrupt 
and  incipient  faults. 


CHAPTER  VI 


Conclusion 


In  [49],  Ljung,  one  of  the  leaders  in  the  field,  argues  that  while  system  identifica¬ 
tion  is  a  mature  area,  there  are  several  important  problems  that  are  not  sufficiently 
understood.  Two  issues  that  he  feels  require  more  investigation  are  parameter  set 
estimation  and  the  tracking  of  time- varying  properties  of  systems  and  signals.  In  this 
dissertation,  we  extended  the  OVE  algorithm  for  parameter  set  estimation  of  linear 
time-invariant  systems  to  allow  for  tracking  of  time-varying  parameters.  Our  expe¬ 
rience  with  the  OVETV  algorithm  challenged  us  to  extend  and  improve  some  of  the 
properties  of  this  algorithm;  several  of  these  modifications  and  adaptations  are  also 
included  in  this  dissertation.  Finally,  we  demonstrated  how  the  OVE  and  OVETV 
algorithms  could  be  utilized  for  fault  detection  and  isolation. 


6.1  Summary 

The  main  result  of  this  dissertation  can  be  found  in  Chapter  II.  In  this  chapter, 
we  combined  the  optimum  volume  time  update  equations  of  [15]  with  the  optimum 
volume  measurement  update  equations  of  [1]  to  develop  the  Optimal  Volume  Ellipsoid 
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algorithm  for  Time- Varying  systems  (OVETV).  Given  a  ]?non  assumptions  on  bounds 
of  the  noise  entering  the  system  and  the  possible  deviation  of  the  parameters  during 
each  time  step,  this  algorithm  is  guaranteed  to  contain  the  true  parameter  as  it  varies 
with  time. 

Building  on  the  development  of  the  OVETV,  we  created  two  algorithms  which 
require  fewer  computations  than  the  optimal  time  update  equations,  thereby  mak¬ 
ing  real-time  implementation  more  feasible  for  many  applications.  These  algorithms 
reduced  the  computational  complexity  of  the  time  update  equations  by  constrain¬ 
ing  the  new  ellipsoid  to  be  parameterized  by  the  previous  ellipsoid  and  a  single  new 
parameter.  One  of  these  algorithms  uses  the  same  parameterization  as  the  scalar 
bound  inflation  method  of  [13].  However,  the  algorithms  developed  in  this  chapter 
are  optimized  for  volume  and  combined  with  the  optimal  volume  measurement  update 
equations  of  [1]. 

In  Chapter  III,  we  applied  the  OVETV  algorithm  to  three  examples.  The  first 
example,  a  simple  first-order  system,  illustrated  the  key  features  of  the  OVETV 
algorithm  and  demonstrated  the  applicability  of  two  alternative  algorithms,  scalar 
addition  and  scalar  multiplication. 

In  the  second  example,  the  effects  of  sampling  and  quantization  on  a  linear  time- 
varying  system,  a  linear  time-varying  circuit,  were  explored.  Theorem  5  was  given 
for  transforming  a  discrete-time  linear  time-varying  state  space  equation  to  a  time- 
varying  difference  equation.  Interestingly  enough,  while  the  goal  of  the  OVETV 
algorithm  is  to  identify  a  parameter  set,  the  center  estimate  of  the  OVETV  algorithm 
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was  able  to  track  the  “true”  parameter  vector  better  than  the  “best”  estimate  using 
WRLS. 

In  the  final  example,  we  applied  the  OVETV  algorithm  to  actual  data  obtained 
from  a  distillation  column.  With  the  “true”  parameter  vector  effectively  unknown,  we 
were  not  able  to  compare  it  with  the  parameter  sets  and  center  estimates.  However, 
the  predicted  output  based  on  the  center  estimate  matched  the  actual  output  very 
closely. 

In  Chapter  IV,  we  discussed  several  algorithm  modifications  and  extensions  to 
the  OVE  and  OVETV  algorithms.  First,  we  showed  how  MISO  systems  could  be 
placed  in  the  OVE  or  OVETV  framework.  Theorem  5  was  extended  from  LTV  SISO 
systems  to  LTV  MISO  systems;  an  extension  to  MIMO  systems  was  also  discussed. 
Results  were  also  presented  showing  the  application  of  the  OVETV  algorithm  within 
the  MISO  framework  to  the  distillation  column  of  Chapter  III. 

Next,  we  showed  how  dependencies  in  parameter  variations  could  be  exploited  to 
reduce  the  number  of  calculations  necessary  to  implement  the  optimal  time  update 
equations  of  Section  2.4.2.  When  these  dependencies  are  present,  the  sizes  of  the 
generalized  eigenvalue  problem  and  the  governing  polynomial  may  be  significantly 
reduced. 

A  square-root  implementation  of  the  OVE  algorithm  is  developed  which  prevents 
the  ellipsoid  orientation  matrix,  Fk,  from  becoming  non-positive  definite  due  to  finite 
precision  calculations.  A  comparison  with  the  direct  implementation  of  the  OVE 
algorithm  showed  that  while  the  square-root  implementation  requires  more  floating 
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point  operations,  the  number  of  operations  is  of  the  same  order  for  each  implementa¬ 
tion.  Furthermore,  in  most  cases  this  increase  in  number  of  computations  will  be  less 
important  than  the  increase  in  numerical  stability  which  is  gained  by  the  square-root 
implementation.  For  the  time- varying  case  it  was  shown  how  the  scalar  addition  and 
scalar  multiplication  algorithms  of  Chapter  II  can  also  be  implemented  in  a  square- 
root  framework. 

The  last  subject  investigated  in  Chapter  IV  was  the  OVE-ISP  input  synthesis 
procedure  of  [21]  and  [11].  First,  the  OVE-ISP  algorithm  was  extended  to  handle 
systems  which  contain  a  known  transportation  lag.  Next,  we  demonstrated  how  OVE- 
ISP  and  OVETV  could  be  applied  to  linear  time-varying  systems.  The  simulation 
results  demonstrated  a  “stabilizing”  property  of  OVE-ISP  procedure  which  is  not 
possessed  by  a  random  input  sequence  and  an  alternative  synthesis  procedure  found 
in  [40].  This  “stabilizing”  property  makes  OVE-ISP  the  preferred  input  scheme  for 
unstable  linear  time-invariant  and  linear  time- varying  systems. 

In  Chapter  V,  we  demonstrated  how  the  OVE  and  OVETV  algorithms  could  be 
used  for  FDI.  Two  different  methods  for  detection  of  failures  were  discussed.  The 
first  method  used  the  consistency  check  which  is  integral  to  the  OVE  and  OVETV 
algorithms.  Because  the  check  is  integral  to  the  OVE  algorithms,  this  detection 
strategy  was  easy  to  implement.  It  was  particularly  adept  at  detecting  abrupt  failures 
in  nominally  slowly-varying  systems. 

An  ellipsoid  intersection  test  was  also  used  to  detect  when  an  ellipsoid  produced 
under  nominal  conditions  fails  to  intersect  with  the  current  ellipsoid  of  the  OVETV 
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algorithm.  The  advantage  of  this  algorithm  lies  in  the  fact  that  it  is  able  to  track  in¬ 
cipient  faults  without  resorting  to  recovery  strategies.  An  alternative  development  of 
the  ellipsoid  intersection  test  equations  of  [46]  was  realized  by  framing  the  intersection 
test  as  a  constrained  optimization  problem. 

Also  proposed  were  two  methods  for  recovering  the  true  parameter  once  a  fault 
causes  an  inconsistency  in  the  OVE  or  OVETV  algorithm.  The  ellipsoid  resetting 
algorithm  resets  the  current  ellipsoid  to  one  large  enough  to  capture  the  true  pa¬ 
rameters.  If  the  system  satisfies  the  nominal  conditions  after  the  ellipsoid  is  reset, 
the  algorithm  is  guaranteed  to  contain  the  true  parameters  until  another  fault  oc¬ 
curs.  Several  different  strategies  were  considered  for  choosing  the  resetting  ellipsoid, 
Ereset-  It  was  shown  that  for  a  particular  form  of  Ereset^  that  as  the  size  of  Ereset  is 
increased,  the  center  estimate  approaches  the  center  estimate  of  our  other  recovery 
strategy,  ellipsoid  projection. 

The  ellipsoid  projection  strategy  results  in  smaller  ellipsoid  volumes  than  the 
ellipsoid  resetting  strategy.  It  combines  the  information  contained  in  the  new  mea¬ 
surements  and  assumes  correct  information  from  the  old  ellipsoid  whenever  ambiguity 
arises.  Because  this  is  not  always  the  case,  repeated  applications  of  the  projection 
algorithm  may  be  required  to  recapture  the  true  parameter.  To  take  advantage  of 
the  guaranteed  containment  properties  of  the  ellipsoid  resetting  algorithm  and  the 
smaller  volumes  of  the  ellipsoid  projection  algorithm,  an  integrated  approach  was 
also  proposed. 


Finally,  during  the  course  of  this  research,  investigation  into  time-varying  linear 
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systems  was  imperative.  In  so  doing,  we  derived  some  very  useful  routines  using 
symbolic  computations  for  linear  time-varying  (LTV)  systems  which  are  given  in 
Appendix  B. 

In  short,  the  major  contributions  of  this  dissertation  are: 

•  development  of  the  OVE  algorithm  for  time- varying  systems, 

•  development  of  alternative  algorithms  for  PSE  of  time-varying  systems,  scalar 
addition  and  scalar  multiplication, 

•  extension  of  OVE  and  OVETV  algorithms  to  MISO  systems, 

•  exploitation  of  dependencies  in  parameter  variations  to  reduce  computations  in 
optimal  time  update  equations, 

•  modification  of  the  OVE  algorithm  using  square-root  approach  for  improved 
numerical  stability, 

•  development  of  PSE  estimation  schemes  for  fault  detection,  and 

•  development  of  schemes  for  parameter  recovery  after  fault  detection. 

6.2  Future  Work 

While  there  are  many  directions  towards  which  this  research  could  turn,  there  are  a 
few  areas  which  seem  particularly  promising.  In  Chapter  IV,  the  simulation  results 
presented  demonstrate  a  “stabilizing”  property  of  OVE-ISP  which  is  not  possessed  by 
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other  input  sequences.  This  stabilizing  property  could  be  potentially  very  important 
for  the  identification  of  unstable  systems  in  the  closed  loop.  However,  more  analysis 
is  needed  to  determine  when  and  if  this  property  can  be  guaranteed. 

At  the  end  of  Chapter  V,  we  discussed  how  the  ellipsoid  intersection  test  could  be 
used  in  conjunction  with  recovery  strategies  for  fault  isolation.  Certainly,  more  work 
is  required  to  explore  the  use  of  PSE  for  fault  isolation  and/or  reconfigurable  control. 
Along  these  lines,  it  would  be  worthwhile  to  explore  how  the  OVE-ISP  procedure 
could  be  modified  to  aid  in  the  isolation  process. 

Most  analytical  redundancy  approaches  for  EDI  involve  either  parameter  estima¬ 
tion  or  state  estimation  schemes.  In  Chapter  V,  we  explored  the  use  of  PSE  for  EDI. 
An  alternative  approach  would  be  to  apply  the  closely  related  set-membership  state 
estimation  schemes  to  the  EDI  problem.  Eor  some  possible  starting  points,  see  [50]. 

Another  area  for  future  investigation  involving  PSE  is  to  use  the  divided-difference 
or  d-operator  instead  of  the  traditional  shift  operator  where  S  =  ^  and  A  is  a 
positive  number.  Traditional  digital  signal  processing  and  control  algorithms  using 
the  shift  operator  often  suffer  from  ill-conditioning  at  high  sampling  rates.  The  ^- 
operator  has  been  shown  to  eliminate  many  of  these  numerical  difficulties  [51,  52].  It 
has  also  allowed  for  a  unification  of  continuous  and  discrete  time  algorithms  [53,  51]. 
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The  following  equations  define  the  parameters  pi  through  ps  used  in  lemma  3. 

Pi  =  -8  -  24  blCal  -  24  -  8  Ca\  -  b\C  -  76  b^^al 

+282  btC^4  -  76  C^albl  -  -  8  b^aj  -  24  6f -  24  a\b{ 

—%b\a\  (A.l) 

P2  =  -17U/CV-1716iV«i"-333fe/CV  +  6  6iW-333  6iV«i" 
-39  +  15  +  6  +  15  -  39  6i®C^ 

+6  6i' V  +  6  +  2  67  +  2  6i®fli®  +  2  C®ai®  (A.2) 

Ps  =  ~  9C\/^V^6i®  -  3C\/^V^6i^ai^ 


+3C"^^\/36i^  +  26i'2 

(A.3) 

P4  =  66i2^ 

P2  +  Ps 

432  6i® 

(A.4) 

Ps  =  66i^  ^ 

P2  -ps 

432  6i® 

(A.5) 
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\Jp4  +P5-  fc/  -  +  5  +  5  ai^i^ 

(A.6) 

6i\/6 

P7  = 

-2Cal  +  +  (^  +  -  2alal  -  2C^a^ 

(A.7) 

P8  = 

:  +  G^Gj  +  a^bl  —  Oi  —  C^b^  +  (^al 

(A.8) 

Appendix  B 


Symbolic  Computations  for  Linear  Time- Varying 

Systems 


B.l  Overview 

During  this  work,  it  has  become  clear  that  the  state  transition  matrix  is  integral  for 
both  analysis  and  control  system  design  of  linear  time-varying  (LTV)  systems.  Un¬ 
fortunately,  a  closed  form  solution  for  the  state  transition  matrix  exists  only  when 
the  LTV  system  satisfies  certain  properties.  In  this  appendix  we  show  how  Maple 
V,  a  computer  algebra  system,  can  be  used  to  calculate  the  state  transition  matrix 
for  several  classes  of  LTV  systems.  By  working  with  symbolic  rather  than  numerical 
data,  computer  algebra  systems  offer  several  advantages  including  the  greater  accu¬ 
racy  achieved  by  the  absence  of  finite  precision  arithmetic  and  the  additional  insight 
obtained  by  maintaining  the  mathematical  information  and  structure.  Examples  and 
applications  to  control  system  design  are  discussed. 
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B.2  Motivation 

Although  the  laws  of  physics  do  not  change  with  time,  time- varying  models  arise  due 
to  special  circumstances  in  the  physical  plant  or  due  to  the  particular  formulation 
of  the  model  [30].  Applications  of  linear  time- varying  (LTV)  systems  include  rocket 
dynamics,  time-varying  linear  circuits,  satellite  systems,  and  pneumatic  actuators. 
The  LTV  structure  is  also  often  assumed  in  adaptive  and  standard  gain- scheduled 
control  systems. 

In  this  appendix,  we  are  interested  in  LTV  systems  of  the  form 

x{t)  =  A{t)x{t) B{i)u(t),  (B.l) 

y(t)  =  C(t)x{t)  + D(t)u{t), 

where  x{t)  €  3^?"  is  the  state  vector,  u{t)  G  3?”^  is  the  control  input,  y{t)  G  3?^  is 
the  system  output,  and  A(it),  B{t),  C'(i),  D{t)  are  matrices  of  appropriate  dimension, 
each  with  (possibly)  time-varying  entries.  The  state  transition  matrix  is  the  unique 
solution  to 

^{t,to)  =  A{t)^{t,to),  $(io,to)  =  A  (B.2) 

where  I  is  the  identity  matrix.  The  state  transition  matrix  is  essential  in  determining 
the  complete  solution,  stability,  controllability,  and  observability  of  (B.l).  It  is  also 
useful  in  design  of  controllers  and  observers  for  (B.l).  Unfortunately,  a  closed  form 
solution  to  (B.2)  exists  only  when  A{t)  satisfies  certain  properties. 

By  working  with  symbolic  rather  than  numerical  data,  computer  algebra  systems 
offer  several  advantages  for  the  controls  engineer,  particularly  for  problems  such  as 
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that  outlined  above.  These  advantages  include  the  greater  accuracy  achieved  by  the 
absence  of  finite  precision  arithmetic  and  the  additional  insight  obtained  by  maintain¬ 
ing  the  mathematical  information  and  structure  [54].  Computer  algebra  systems  also 
handle  complex  mathematical  manipulations,  which  if  done  by  hand  would  be  tedious 
and  error  prone  [55].  One  particular  computer  algebra  system,  which  we  use  in  this 
work,  is  Maple  V,  which  is  a  powerful  software  package  for  symbolic  and  numeric 
computation.  Maple  V  includes  a  large  library  of  functions,  programming  capability, 
interactive  graphics,  a  worksheet  interface,  and  an  online  help  facility.  It  is  available 
on  many  computer  platforms  including  MS-DOS,  Windows,  Macintosh,  NeXT,  DEC, 
Sun  and  other  UNIX  workstations. 

Maple  V  can  easily  find  the  solution  of  (B.2)  when  A{t)  is  constant  or  satisfies 
a  well  known  commutative  property.  This  appendix  describes  routines  written  by 
the  authors  to  expand  the  class  of  systems  for  which  Maple  V  can  easily  calculate 
^{t,to)  to  include  systems  where  A{t)  is  triangular  or  A{t)  satisfies  certain  bracket 
properties.  Routines  have  also  been  written  to  calculate  the  general  solution  of  (B.2) 
to  any  arbitrary  order  and  to  check  whether  a  given  satisfies  (B.2). 


B.3  State  Transition  Matrix  Properties 

The  state  transition  matrix  is  an  integral  component  in  the  study  of  LTV  systems  of 
the  form  given  by  (B.l).  It  is  used  for  determining  the  complete  solution,  stability, 
controllability,  and  observability.  It  can  also  be  used  in  the  design  of  controllers  and 
observers  for  (B.l).  In  this  section  we  will  discuss  these  uses  along  with  some  of  the 
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properties  of  the  state  transition  matrix.  Since  this  material  is  somewhat  standard  in 
the  linear  systems  literature,  we  offer  here  only  a  brief  overview,  and  present  relevant 
theorems  in  the  Addendum. 

The  state  transition  matrix,  $(t,to)5  satisfies 

i(t,to)  =  A(t)$(t,to),  (B.3) 

and  has  the  following  important  properties  [56]: 

^(t,t)  =  I  (B.4) 

(B.5) 

^(t2,to)  =  ^t2,ti)^{ti,to).  (B.6) 

Stability  of  the  homogeneous  system, 

x{t)  =  A{t)x{t),  (B.7) 

whose  solution  is  given  by 

x(t)  =  $(t,to)xo,  (B.8) 

where  xq  =  ^(^0),  can  be  determined  from  the  state  transition  matrix,  according  to 
well  known  stability  theorems  [56]  (see  Addendum).  The  necessary  and  sufficient 
conditions  on  $(t,  to)  for  stability  are  summarized  in  Table  7. 

It  is  easy  to  verify  that  the  solution  to  the  non- homogeneous  system  (B.l)  is  given 
by 

x{t)  =  ^{t,to)xo+  f  $(t,T)B(T)u(r)dT 

Jto 

y{t)  =  C(t)^{t,to)xo  +  C{t)  f  ^{t,T)B{T)u{T)dT  +  D{t)u{t).  (B.9) 

Jto 
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Table  7:  Stability  bounds  on  $(<,  to) 


Stability  Result 

Necessary  and  Sufficient  Condition 

Stable  in  the  sense  of  Lyapunov  at  to 

l|$(Ufo)||  <  k{io)  <  oo 

Uniformly  stable  in  the  sense  of  Lyapunov 

ii$(f,<o)il  <  k  <  oo 

Asymptotically  stable  at  to 

ii$(t,fo)ii  <  k{to)  <  oo 
ii$(Uto)ii  — >  0  as  f  — »■  oo 

Uniformly  asymptotically  stable 

mt,to)\\  < 

To  guarantee  that  the  system  can  be  driven  from  one  state  xq  to  another  state 
Xi  with  an  input  u(t),  it  is  necessary  to  show  that  the  system  is  controllable.  The 
LTV  system  (B.l)  is  said  to  be  controllable  if  given  any  xq  there  exists  an  input 
u(t)[(o,4i]  such  that  x{ti)  =  0.  Controllability  of  (B.l)  can  be  determined  from  the 
state  transition  matrix  according  to  a  well  known  theorem  [30]  (see  Addendum). 

To  guarantee  that  the  system  state  x{t)  can  be  estimated  from  the  system  output 
y{t),  it  it  necessary  to  show  that  the  system  is  observable.  The  LTV  system  (B.l) 
is  said  to  be  observable  on  [to)^i]  if  the  initial  state  Xo  is  uniquely  determined  by 
the  output  y{t)  for  t  e  Observability  of  (B.l)  can  be  determined  from  the 

state  transition  matrix  according  to  a  well  known  theorem  [30]  (see  Addendum).  The 
controllability  and  observability  grammians,  W {to,  h)  and  M{to,  h),  respectively,  (see 
Addendum)  can  also  be  used  in  the  design  of  controllers  and  observers  for  (B.l). 

It  is  clear  that  the  state  transition  matrix  is  important  for  studying  stability, 
controllability,  and  observability  of  (B.l).  Calculation  of  the  state  transition  matrix 
for  linear,  time-invariant  systems  is  a  straight  forward  task.  Unfortunately,  for  lin- 
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ear  time-varying  systems,  it  is  often  difficult  if  not  impossible  to  calculate  the  state 
transition  matrix. 

B.4  Calculating  the  State  Transition  Matrix 

In  general,  a  closed  form  solution  for  does  not  exist.  In  this  section,  several 

classes  of  systems  for  which  can  be  calculated  in  closed  form  are  discussed. 

Maple  V  routines  which  aid  in  this  calculation  are  discussed.  After  summarizing 
classes  for  which  the  state  transition  can  be  calculated,  two  decomposition  schemes 
will  be  discussed  which  expand  on  these  classes. 

Before  examining  these  classes,  some  comments  about  Maple  V  should  be  given. 
As  mentioned  earlier.  Maple  V  contains  a  large  library  of  functions.  While  many  of 
these  functions  are  internal  to  Maple  V  and  available  immediately,  other  functions  are 
grouped  together  as  packages  which  must  be  loaded  into  Maple  V  before  using  those 
commands.  One  such  package,  the  linear  algebra  package,  contains  many  common 
functions  for  working  with  vectors  and  matrices  [57].  Several  of  the  functions  that  we 
use  in  this  work  are  given  in  Table  8.  While  many  of  Maple  V’s  other  commands  are 
self-explanatory,  the  procedure  map  also  requires  some  discussion.  The  procedure 
map(f,  A,  arg2,  argS,  ...,  argn)  is  used  to  apply  a  function,  /,  with  multiple 
arguments  to  each  component  of  an  expression  (or  matrix),  A. 

In  addition  to  the  linear  algebra  package,  it  is  assumed  that  several  procedures, 
which  were  written  by  the  authors  and  listed  in  Table  9,  are  also  loaded.  The  eye 
and  zeros  procedures  are  modeled  after  the  Matlab  functions  of  the  same  names  and 
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Table  8:  Maple  V  linear  algebra  procedures 


linsolve(A,  b) 

Solves  Ax  =  b 

evalm(-) 

Evaluates  matrix  expression 

add(A,.B) 

Adds  matrices  A  and  B 

multiply(A,  B) 

Multiplies  matrices  A  and  B 

exponential(A,  t) 

Computes  matrix  exponential, 

transpose(A) 

Matrix  transpose,  A^ 

eigenvals(A) 

Computes  eigenvalues  of  matrix  A 

return  the  identity  matrix  and  a  matrix  of  zeros,  respectively.  The  msubs  procedure 
substitutes  expressions  into  matrices.  The  Kronecker  product  of  A  =  [ajj]  E 
and  B  E  is  defined  to  be 


A®B 


aiiB 


O’lnB 


'  ’  ‘  (^mnB 

where  A®  B  E  The  vector  vec  A  is  defined  to  be 


(B.IO) 


vec  A  = 


Ol 


(B.ll) 


L  «n  J 

where  A  E  vec  A  E  and  the  a,-  are  the  columns  of  the  matrix  A.  The 

command  invec(vec(A)  ,m,  n)  will  return  the  matrix  A.  The  procedures  kron  and 
vec  transform  matrix  equations  of  the  form 


AXB  =  C 


(B.12) 


where  C  €  and  X  G  is  unknown,  to  the  equivalent  system  of  equations 


®  A)vec  X  —  vec  C. 


(B.13) 
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Table  9:  Maple  V  utility  procedures 


eye(-) 

Generates  identity  matrix 

zeros(-) 

Generates  matrix  of  zeros 

msubs(y4,  si, ...,  Sn) 

Substitutes  expressions  into  a  matrix 

kron(A,B) 

Kronecker  product  A®  B 

vec(A) 

vec  A 

invec(a,m,n) 

inverse  of  vec  A 

This  equation  can  be  solved  for  vec  X  by  linsolve(D,  vec(C))  where  D  =  ®  A 

[58]. 

With  these  tools,  we  can  now  study  the  solution  of  the  general  problem  (B.2). 
While  in  general  there  is  no  closed  form  solution  to  (B.2),  the  solution  can  be  expressed 
in  terms  of  the  Peano-Baker  series  [30] 

$(t,to)  =  7+/  A(ri)dri  +  /  A{ti)  [  A(r2)dr2dri  +  •  •  • .  (B-14) 

JtQ  v  to  *'tQ 

The  following  procedure,  peano(A,  n,  tO),  allows  us  to  calculate  this  series  to  any 
prespecified  order  n: 


peano : = 
proc(A,n,tO) 
local  wl,w2,i,tp; 

if  not  type(n,posint)  then  ERR0R('n  must  be  positive  integer') 
elif  not  type (A, matrix)  then  ERR0R('A  must  be  matrix') 
elif  not  type (to, scalar)  then  ERROR('tO  must  be  a  scalar') 

fi; 

w2  ;=  eye (A) ; 
wl  ;=  evalm(w2) ; 
for  i  to  n-1  do 
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w2  :=  mapCint .multiply (A, w2) ,t  =  tO  . .  tp) ; 
w2  :=  msubs(msubs(w2,t  =  t.i),tp  =  t) ; 
wl  :=  add(wl,w2) 

od; 

RETURN (") 

end 


For  time-invariant  systems,  that  is,  when  A{t)  =  A  is  a  constant  matrix,  the  state 
transition  matrix  is  given  as 

When  A  is  time-invariant.  Maple  V  can  easily  calculate  ^{t,to)  using  the  Maple  V 
function  exponential(A,  t-tO). 

If  A{t)  satisfies  the  commutative  property, 

A{t)  A{h)dU^  =  (£  A{h)dh^  A{t),  (B.16) 

for  all  t,  to,  the  state  transition  matrix  is  given  as 

=  (B.17) 

Again,  Maple  V  can  easily  calculate  using  the  command  exponential  (map 

(int,msubs(A,t=tl),  tl=t0..t)). 

Checking  (B.16)  is  equivalent  to  checking  if 

A{h)A{t2)  =  A{h)A{h),  Vti,  h  (B.18) 


is  satisfied  [59].  The  commutative  condition  can  also  be  checked  by  decomposing  A{t) 
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into  the  form 

A{t)  =  j2c^i{t)Ai  (B.19) 

i=i 

where  the  ai(t)’s  are  linearly  independent  time  functions  and  the  Aj’s  are  constant 
matrices,  and  then  checking  if  the  A,  ’s  form  a  commuting  family,  that  is  AiAj  =  AjAi, 
Vi,j.  When  the  system  is  commutative,  the  state  transition  matrix  is  given  by 

(B.20) 

i=l 

where  0i{t,to)  =  J^^ai{T)dT  [60].  The  solution  of  (B.18)  and  (B.20)  are  equivalent; 
however,  the  solution  of  (B.20)  may  result  in  a  simpler  form.  For  another  decomposi¬ 
tion  method,  see  [61].  Clearly,  constant  matrices  and  diagonal  matrices  both  satisfy 
the  commutative  property. 

If  A{t)  exists  and  there  is  a  constant  matrix  Ai  which  satisfies 

^A{t)  =  AiA{t)  -  A{t)Ai  Vt  >  to,  (B.21) 

then 

$(t,to)  =  Vt  >  to,  (B.22) 

where  A2  =  .di(to)  —  Ai  [62].  If  A{t)  satisfies  (B.21),  it  is  said  to  be  a  member  of  the 
Ai  class.  A  necessary  condition  to  satisfy  (B.21)  is  that  the  eigenvalues  of  A(t)  are 
time-invariant.  Stability  for  members  of  Ai  can  be  completely  determined  from  the 
eigenvalues  of  Ai  and  A2. 
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The  Kronecker  product  is  used  to  transform  (B.21)  to 

{A^  ®  I  —  I  ®  A)vec  Ai  =  vec  A,  (B.23) 

which  can  be  solved  for  vec  Ai  using  the  Maple  V  routine  linsolve.  The  following 
procedure,  fal(A,  AO,  AI,  A2),  transforms  (B.21)  to  (B.23)  and  returns  Ao,  Ai, 
and  A2: 


fal:  = 

proc ( A , AO , AI , A2) 
local  ad,Ad,al ,A1 ,n; 

n  :=  linalgCrowdim] (A) ; 

Ad  :=  mapCdiff ,A,t) ; 
ad  :=  vec (Ad); 

AI  :=  add(kron(linalg[traiispose]  (A)  ,eye(A))  ,-kron(eye(A)  ,A))  ; 
al  : =  linalg [linsolve] (AI , ad) ; 

if  al  =  NULL  then  ERR0R('no  solution  possible')  fi; 

AO  :=  msubs(A,t  =  0); 

Al  :=  invec(al,n,n) ; 

A2  :=  evalm(AO-Al) 

end 


We  note  that  multiple  solutions  are  possible;  the  user  should,  if  possible,  choose  the 
free  parameters  such  that  Ai  and  A2  are  time-invariant.  If  a  time-invariant  solution 
exists,  ^{t,to)  is  given  by  multiply  (exponential  (Al,  t-tO),  exponential  A2,  t- 


to)). 

The  class  Ai  can  be  extended  as  follows  [63].  If  A{t)  and  h{t)  exist,  where  h{t)  is 
a  nonzero  scalar  time  function,  and  there  is  a  constant  matrix  Ai  which  satisfies 


m 

h{t) 


d 

dt 


=  A-iA{t)  —  A{t)Ai  yt  >  to, 


(B.24) 
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then 

to)  =  >  to,  (B.25) 

where  g{t)  =  /i(r)c?T  and  A2  =  -  Ai.  If  A{t)  satisfies  (B.24),  it  is 

said  to  be  a  member  of  the  Ah  class.  A  necessary  condition  to  satisfy  (B.24)  is 
that  the  eigenvalues  of  A{t)  are  scalar  multiples  of  h{t).  This  can  be  checked  using 
eigenvals(A).  When  h{t)  =  1  this  class  reverts  to  the  Ai  class,  and  when  h{t)  =  1/t 
the  Euler-type  differential  equation  can  be  solved  [64]. 

The  following  procedure,  fah(A,  AhO,  AI,  A2,  tO,  h,  g),  transforms  (B.24) 
using  the  Kronecker  product  and  returns  Ai,  A2  and  g: 


f  ah:  = 

proc(A,Ah.0,Al,A2,t0,h,g) 
local  ad,Ad,al,Al,n,Ah,tl; 

g  :=  eval(int(subs(t  =  tl,h),tl  =  tO  . .  t)); 
n  :=  linalgCrowdim] (A) ; 

Ah  :=  evalm(l/h*A) ; 

Ad  :=  map(diff ,Ah,t) ; 
ad  :=  vec(Ad); 

AI  :=  add(kron(linalg [transpose] (A) ,eye(A)) ,-kron(eye(A) , A)) ; 
al  :=  linalgClinsolve] (Al,ad) ; 

if  al  =  NULL  then  ERR0R(‘no  solution  possible')  fi; 

AhO  :=  map(limit ,Ah,t  =  tO) ; 

Al  :=  invec(al,n,n) ; 

A2  :=  evalm(AhO-Al) 

end 


Again,  if  a  time-invariant  solution  for  Ai  and  A2  exists,  $(t,to)  is  given  by  multiply( 
exponential(Al,  g),  exponential(A2,  g)). 
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Often  situations  arise  where  A{t)  is  triangular.  If  A(t)  =  [aij(t)]  is  lower  triangular, 
that  is  aij{t)  =  0,  Vj  >  i,  then  the  state  transition  matrix  has  elements  given  by 


t  <  j 

i=j 


Jtl  WdT 
^  to)dT  i  >  j, 


(B.26) 


i— 1 


k=j 


where  The  following  procedure,  Itriangle(A),  calculates 

for  lower  triangular  matrices: 


ltriangle:= 

proc(A) 

local  phi,n,i,sumk,k,j ,1; 
n  :=  rowdimCA); 
phi  :=  map(0,A) ; 
for  i  to  n  do 

phi[i,i]  :=  exp ( int( subs (t  =  tl,A[i,i] )  ,tl  =  tO  . .  t)); 
for  i  to  i-1  do 
1  :=  i-j+1; 

sumk  :=  subs(t  =  t .l,sum(A[i,k]*phi [k, j] ,k  =  j  ..  i-1)); 
phi[i,j]  :=  phi[i,i]* 

int(subs(t  =  t0,subs(t0  =  t  .l,phi[i,i]))*suink,t.l  =  t0..t) 
od 

od; 

map (simplify, phi) ; 

RETURNC") 

end 


If  A{t)  is  upper  triangular,  that  is  aij{t)  =  0,  Vz  >  j,  the  state  transition  has 
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elements  given  by 


to) 


r‘  aii{T)dT 

e  ‘0 


t  >  J 

i^j 


(B.27) 


/  4>n{tyT)  Clik{T)(f>kj{T,to)dT  i  <  j. 
•^*0  k=i+l 


Note  that  a  diagonal  matrix  is  both  lower  and  upper  triangular.  The  following  pro¬ 
cedure,  utriangle(A),  calculates  $(t,to)  for  upper  triangular  matrices; 


utriangle  :  = 
proc(A) 

local  phi,n,i,suink,k,j  ,1; 
n  :=  rowdim(A); 
phi  :=  map(0,A) ; 
for  i  from  n  by  -1  to  1  do 

phi[i,i]  :=  exp (int (subs (t  =  tl,A[i,i]),tl  =  tO  ..  t)); 
for  1  from  i+1  to  n  do 
1  :=  j-i+1; 

sumk  ;=  subs(t  =  t .l,sum(A[i,k]*phi [k, j] ,k  =  i+1  ..  j)); 
phi[i,j]  :=  phi[i,i]* 

int(subs(t  =  t0,subs(t0  =  t .l,phi [i , i] ))*sumk,t . 1  =  t0..t) 
od 
od; 

map (simplify, phi) ; 

RETURNC) 

end 


This  concept  can  be  extended  to  block  triangular  matrices.  Let  A{t)  =  [A,y(t)], 
where  each  Aij{t)  has  dimensions  n,-  x  ny  and  #(t,to)  =  where 

has  dimensions  n,  x  ny  If  A{t)  is  lower  block  triangular,  that  is  Aij{t)  =  0,  Vj  >  i, 
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then  the  state  transition  matrix  has  blocks  given  by 


^ij{t,to)  =  \ 


$ 


^  <  j 


t=j 


(B.28) 


pt  *  1 

/  ^ii{t,T)y]Aik{T)^kj{T,to)dT  i>j. 
•^‘0  k=j 


Likewise,  if  A{t)  is  upper  block  triangular,  that  is  Aij{t)  =  0,  Vi  >  j,  then  the  state 
transition  has  blocks  given  by 


to)  — 


0  i>  j 

$n(t)  to)  i  j 

/  $ii(t,T)  Y,  Aik{T)^kj{T,to)dT  i<j. 

•'tn  J _ •  I  1 


(B.29) 


While  (B.28)  and  (B.29)  hold  for  all  lower  and  upper  block  triangular  matrices,  re¬ 
spectively,  #(t,to)  can  be  calculated  explicitly  only  when  the  $ij(t,to)  are  known. 
Therefore,  if  every  Ai(t)  is  from  a  class  of  matrices  from  which  $ji(t,to)  can  be  calcu¬ 
lated,  $(t,  to)  can  also  be  calculated.  Note  that  a  block  diagonal  matrix  is  both  lower 
and  upper  block  triangular. 

Two  decomposition  schemes  exist  which  expand  the  class  of  systems  for  which 
$(t,  to)  can  be  calculated.  The  first  one  can  be  found  in  [65].  Let  ^(t)  be  decomposed 
as 

m 

A(t)  =  ^A,{t) 

i=l 


(B.30) 
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such  that 

=  Fi(t)^i(t,0), 

i  =  l,2,...,m  (B.31) 

H0,0)  =  I, 

is  solvable  where 


Fi(t)  =  Ti_\(t)Ai(t)Ti-i(t)  i  =  1,2,.. .  ,m, 


{B.32) 


and 


Ti{t)  =  ^i{t,0)Ti-i{t)  i  =  l,2,...,m-l, 
Toit)  =  I. 


(B.33) 


Then  we  have  that 

m 

«((,0)  =  n*i((,0).  (B.34) 

»=1 

The  procedures  discussed  to  this  point  are  useful  in  constructing  procedures  for  com¬ 
puting  this  decomposition.  It  is  interesting  to  note  that  this  decomposition  scheme 
can  be  used  to  show  that  any  arbitrary  A{t)  can  be  decomposed  into  two  normal 
systems  [66].  Consequently,  if  closed  form  solutions  exist  for  these  normal  systems,  a 
closed  form  solution  exists  for  any  arbitrary  A(t). 

A  second  decomposition  scheme  can  be  found  in  [67]  and  [68].  Differentiating 
(B.7)  results  in 


X  =  (A(t)  -f  A^{t))x{t). 


(B.35) 


If 


A(t)  -1-  A^{t)  —  B{t)A(t), 


(B.36) 
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for  some  B(t),  then  (B.35)  can  be  rewritten  as 


X  =  B{t)x. 


(B.37) 


If  a  solution  for  B{t)  exists,  the  following  procedure,  hemami(A),  returns  B{t): 


hemami := 
proc(A) 

local  ad,Ad,b,B,n,Al; 

n  :=  linalgCrowdim] (A) ; 

Ad  :=  mapCdiff ,A,t) ; 
ad  :=  vec(evalm(Ad+A''2))  ; 

A1  :=  kron(linalg[transpose] (A) ,eye(A)) ; 
b  ;=  linalgClinsolve] (Al,ad) ; 

if  b  =  NULL  then  ERR0R('no  solution  possible')  fi; 
B  :=  invec(b,n,n) 

end 

If 


$i(0,0)  =  I  (B.38) 

is  solvable  then 

#(i,io)=  ^i{T,to)dT^  A{to)  +  I.  (B.39) 

Assuming  Phil  has  been  found,  is  given  by  add(multiply(map(int,  msubs 

(Phil,  t=tl),  tl=t0..t),  msubs(A,  t=tO)),  eye(A)).  If  (B.38)  is  not  solvable 
using  any  of  the  previous  procedures,  (B.37)  can  be  differentiated.  If  a  solution  to 
the  Riccati  equation  B{t)  +  B^{t)  =  C{t)B{t)  exists,  this  procedure  can  be  repeated 
until  $(t,to)  is  found,  or  until  the  Riccati  equation  has  no  solution. 
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Lastly,  we  discuss  the  procedure  used  to  verify  that  a  given  satisfies  (B.2). 

From  (B.2),  we  know  that 


A{t)  =  $  \t,to)^{t,to)  (B.40) 

^{to,to)  =  I. 

The  following  procedure  check(Phi)  returns  ^~''-{t,to)^{t,to)  and  $(to,to): 


check : = 
proc(phi,t,tO) 
local  a, Id; 

a  :=  multiply(map(diff , phi, t) .inverse (phi)) ; 
Id  :=  map (limit, phi, t  =  tO) ; 

RETURN (map (simplify, Id) , map (simplify, a)) 

end 


In  this  section,  several  classes  for  which  the  state  transition  matrix  can  be  calcu¬ 
lated  were  given.  These  include  constant  matrices,  matrices  which  satisfy  the  commu¬ 
tative  condition,  the  Ai  class,  the  Ah  class,  triangular  matrices,  and  block  triangular 
matrices  whose  block  diagonal  matrices  are  solvable.  Two  decomposition  schemes 
which  extend  these  classes  were  also  given.  Maple  V  procedures  were  given  which  aid 


in  this  calculation. 
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B.5  Examples 

In  this  section,  we  will  apply  the  Maple  V  procedures  discussed  earlier  to  two  LTV 
systems.  In  exposition  of  the  examples,  we  include  diary-like  output  from  the  Maple 
V  implementation  to  assist  the  reader  in  reproducing  these  results. 

Example  1 

Consider  Euler’s  equations  for  the  rigid  body  [30]: 

h 

W2(^)  =  “7 - “^3(tVl(^)  +  (B-41) 

12  12 

^3(0  “  +  ■7“W3(0 

11  13 

where  the  uji  are  the  angular  velocities,  the  /,•  are  the  principal  moments  of  inertia,  and 
the  Ui  are  the  applied  torques.  If  we  let  E  =  I2  =  I  (symmetrical  body)  and  linearize 
about  the  nominal  trajectory,  ui  =  U2  =  «3  =  0,  wi(t)  =  sm{wgt),  W2(t)  =  cos{wgt), 
and  ui3{t)  =  w,  where  g  =  we  obtain  the  state  equation 


6(t)  =  A{t)Q{t)  -b  B{t)u{t\  (B.42) 

where  u{t)  =  vo{t)  -  cb(i).  In  this  case,  the  matrix  A{t)  can  be  found  according  to 

>  A:=matrix([[0,g*w,g*cos(w*g*t)] , [-g*w,0 ,-g*sin(w*g*t)] , 

>  [0,0,0]]); 


[0  gw 
[ 

A  :=  [  -  g  w  0 

[ 

[  0 


g  cos(w  g  t) 

-  g  sin(w  g  t) 


] 

] 

] 

] 

] 


0 


0 
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The  eigenvalues  of  A{t),  given  by 
>  eigenvals(A) ; 

221/2  22  1/2 
0,  (-  g  w  )  ,  -  g  w  ) 

are  time-invariant,  which  suggests  that  A{t)  may  be  in  the  Ai  class.  The  procedure 
call  fal(A,  AO,  Al,  A2,  0)  returns  Ai  and  A2  with  Ai  given  below  as 


>  fal(A,A0,Al,A2,0)  : 

>  print (Al) ; 


[  cos(w  g  t)t4 

[  t5  t4  - + 

[  w 

[ 

[  -sin(w  g  t)t4 

[-t4  t5  - + 

C  w 

C 

[  0  0 


sin(w  g  t)t5  sin(w  g  t)t9  ] 

-  g  cos(w  g  t)] 

w  w  ] 

] 

cos(w  g  t)t5  cos(w  g  t)t9  ] 

-  +  g  sin(w  g  t)] 

w  w  ] 

] 

t9  ] 


where  t4,  t^,  and  tg  are  parameters  to  be  specified.  These  parameters  are  chosen  such 
that  Al  and  A2  are  time-invariant,  according  to 


>  Al :=msubs(Al ,t4=w*g,t5=0,t9=0) ; 

[  0  g  w  0  ] 

[  ] 

Al:=[-gw  0  0] 

[  ] 

[  0  0  0  ] 


>  A2:=msubs(A2,t4=w*g,t5=0,t9=0) ; 
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1 — 1 

0 

0 

1 — 1 

bO 

1 — 1 

] 

A2  :=  [  0 

0 

0 

1 _ 1 

[ 

1 — 1 

1 — 1 

0 

0 

1 — 1 

0 

Now  0)  can  be  calculated  as  in  (B.22)  using 

>  phi ;=map(simplify,multiply(exponential(Al,t) , 

>  exponential(A2,t))) ; 

[cos(wgt)  sin(wgt)  cos(wgt)  t  g] 

[  ] 

phi :=[-sin(wgt)  cos(wgt)  -sin(wgt)  g  t] 

C  ] 

[0  0  1  ] 

Finally,  checking  our  results  (according  to  (B.40)) 


>  check(phi,t,0) ; 


Cl 

0 

0] 

[  0  gw 

g  cos(w  g  t)  ] 

[ 

] 

c 

] 

[0 

1 

0], 

C-  g  w  0 

-  g  sin(w  g  t)] 

[ 

] 

c 

] 

[0 

0 

1] 

[  0  0 

0  ] 

The  identity  matrix  and  the  matrix  A{t)  are  returned  as  expected.  Note  that  by 
examining  $(t,0),  we  see  that  the  nominal  trajectory  is  unstable  (see  Table  7).  The 
elements  ^i3(t,0)  and  (l)23{t,0)  of  $(t,0)  are  clearly  unbounded. 

Example  2 


Again  consider  Euler’s  equations  for  a  symmetrical  body,  that  is  Ii  =  I2  =  I. 
Setting  ui  =  U2  =  0,  the  equations  can  be  written  in  the  form  of  a  quasi-linear 
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parameter  varying  system  [69]: 

0  gz{t)  0  0 

~  0  0  0  ^3(^)  (B.43) 

z(()  J  [  0  0  0  J  [  4*)  J  [  5 . 

where  z[t)  =  cn3{t).  A  linear  parameter  varying  (LPV)  system  depends  on  a  time- 
varying  parameter  rather  than  explicitly  on  time.  This  parameter  is  known  at  time 
t,  but  not  necessarily  a  priori  [69]. 

The  control  structure  U3{t)  =  —l3k{z  —  Zref)  is  proposed  to  drive  z{t)  to  Zrej  and 
maintain  stability  of  the  system  where  fc  is  a  positive  constant.  Substituting  u^it) 
into  (B.43),  we  obtain  the  following  closed  loop  matrix  for  A(t), 

>  A:=matrix([[0,g*z(t) ,0] , [-g*z(t) ,0,0] , [0,0,-k]]) ; 

[  0  g  z(t)  0  ] 

[  ] 

A  :=  C  -  g  z(t)  0  0  ] 

C  ] 

[  0  0  -  k  ] 

Taking  the  Peano-Baker  series  to  the  5th  order  results  in 

>  peano(A,5,0) ; 

2  4  3 

[1-1/2  %1  +  1/24  yj  ,  •/.2  -  1/6  y.2  ,  0] 

3  2  4 

[y.1-1/6  y.i  ,  1-1/2  y.2  +  1/24  y.2  ,  0] 

2  2  3  3  4  4 

1-k  t+1/2  k  t  -l/6k  t  +l/24k  t  ] 


[0, 


0, 
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t 

/ 

I 

y.l  :=  I  "  g  zCtl)  dtl 

I 

/ 

0 

t 

/ 

I 

7.2  :=  I  g  z(tl)  dtl 

I 

/ 

0 
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By  examining  the  Peano-Baker  series  and  checking  a  table  of  Taylor  series,  it  is 
proposed  that  $(t,0)  is  of  the  following  form: 

>  h:=iiit(g*z(tl)  ,tl=0.  .t)  : 

>  phi :=matrix( [[cos(h) ,sin(h) ,0] , [-sin(h) ,cos(h) ,0] , 

>  [0,0,exp(-k*t)]]) ; 

[  cos(y.l)  sin(y,l) 

[ 

phi:  =  [  -  sin(y.l)  cos (7,1) 

[ 

[  0  0 

t 
/ 

I 

y.l  :=  I  g  z(tl)  dtl 

I 

/ 

0 


0  ] 

] 

0  ] 

3 

exp(-  k  t)  ] 


Checking  our  proposed 
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>  check(phi,t ,0) ; 


[  1 

0 

0  ]  [ 

0 

g  z(t) 

0  ] 

[ 

]  [ 

] 

C  0 

1 

0  ],  [ 

-  g  z(t) 

0 

0  ] 

1 — 1 

]  c 

] 

0 

1 — 1 

0 

1  ]  [ 

0 

0 

-  k  ] 

where  the  identity  matrix  and  A{t)  are  returned  as  expected.  By  examining  0), 
it  can  be  seen  that  oji  and  u)2  are  stable  and  that  2;  is  asymptotically  stable  as  desired. 


B.6  Summary 

In  this  appendix,  we  have  examined  how  Maple  V  could  be  used  to  calculate  the 
state  transition  matrix  for  several  classes  of  LTV  systems.  We  motivated  the  problem 
by  arguing  that  the  state  transition  matrix  is  important  for  understanding  stability, 
controllability,  and  observability  of  LTV  systems.  Classes  of  systems  for  which  the 
state  transition  matrix  can  be  calculated  include  constant  matrices,  matrices  which 
satisfy  a  commutative  condition,  the  Ai  class,  the  Ah  class,  triangular  matrices, 
and  block  triangular  matrices  whose  block  diagonal  matrices  are  solvable.  Maple  V 
procedures  to  aid  in  this  calculation  were  given.  Finally,  examples  showing  how  the 
state  transition  matrix  can  be  calculated  and  used  for  analysis  and  control  of  LTV 
systems  were  discussed.  Clearly,  this  is  one  of  many  applications  where  symbolic 
software  packages  can  be  used  to  aid  the  controls  engineer. 
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Addendum 

In  the  following  we  list  several  well  known  Theorems  (involving  the  state  transition 
matrix),  without  proof,  which  are  referenced  in  the  text. 


1.  Stability  [56]: 

Theorem  10  Every  equilibrium  state  of  (B.7)  is  stable  in  the  sense  of  Lya¬ 
punov  at  to  if  and  only  if  there  exists  some  constant  k  which  depends  on  to  such 
that  |l$(t,to)||  <  k  <  oo  for  all  t  >  to-  If  k  is  independent  of  to,  it  is  uniformly 
stable  in  the  sense  of  Lyapunov. 

2.  Asymptotic  Stability  [56]: 

Theorem  11  The  zero  state  of  (B.7)  is  asymptotically  stable  at  to  if  and  only 
if  ||$(t,io)ll  <  K^o)  <  oo  <^nd  ll$(t,to)||  ^  0  as  t  oo.  The  zero  state  is 
uniformly  asymptotically  stable  over  [0,  oo)  if  and  only  if  there  exist  positive 
numbers  ki  and  k^  such  that  ||$(t,<o)||  <  any  to>0  and  for  all 

t>to. 


3.  Controllability  [30]: 
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Theorem  12  The  LTV  system  (B.l)  is  controllable  on  if  and  only  if  the 

controllability  grammian 

W{to,ti)=f  ^{to,t)B{t)B'^{t)^'^{to,t)dt  (B.44) 

J  tQ 

is  invertible. 

4.  Observability  [30]: 

Theorem  13  The  LTV  system  (B.l)  is  observable  on  if  and  only  if  the 

observability  grammian 

M{to,h)=  r  ^'^{to,t)C^{t)C{t)^to,t)dt  (B.45) 

Jto 


is  invertible. 
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